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ABSTRACT 

A  computer  simulation  of  a  correlator  receiver  was 
developed  and  exercised  to  study  the  impact  of  a  model-based 
signal  processing  algorithm  on  the  detection  of  transmitted 
CW  and  LFM  pulse  acoustic  signals  incident  on  a  planar  array 
of  electroacoustic  transducers.   The  model  of  the  ocean 
communication  channel  incorporates  a  space-variant  sound 
speed  profile.   The  transducer  output  electrical  signals  are 
cophased   by  an  FFT  beamformer  via  phase  weighting,  and 
summed  to  form  a  total  array  output  signal.   The  total  array 
output  signal  is  correlated  with  a  delayed  replica  of  the 
transmit  waveform  and  compared  to  a  Neyman-Pearson  thres- 
hold.  Receiver  performance  is  measured  using  a  Monte  Carlo 
technique  to  estimate  the  probability  of  detection  for  a 
fixed  probability  of  false  alarm  versus  the  signal-to-noise 
ratio  at  the  input  of  a  single  transducer.   White,  zero-mean, 
Gaussian  transducer  noise  is  assumed  to  facilitate  compari- 
son between  theoretical  and  simulated  performance.   Results 
indicate  that  model-based  signal  processing  provides  signi- 
ficant improvement  of  receiver  performance. 
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I.   INTRODUCTION 

Model-based  signal  processing  is  described  by  Mendel 
[Ref.  1]  as  an  approach  that  exploits  knowledge  of  the 
underlying  physics  of  a  problem  to  develop  signal  processing 
algorithms.   Use  of  the  approach  implies  that  some  a  priori 
knowledge  exists  regarding  the  problem  under  consideration. 
In  the  case  of  an  underwater  acoustic  communication  problem, 
such  a  model  has  been  developed  by  Ziomek  [Ref.  2]. 

Ziomek  derived  a  time-invariant,  space-variant  transfer 
function  of  the  ocean  volume  based  on  the  WKB  approximation, 
which  is  an  approximate  solution  of  the  linear,  inhomogeneous , 
scalar  wave  equation  describing  the  propagation  of  small- 
amplitude  acoustic  pressure  waves  when  the  speed  of  sound  is 
a  function  of  depth.   Based  on  the  transfer  function  of  the 
ocean  volume,  Ziomek  [Ref.  2:pp.  257-261]  also  derived  an 
equation  describing  the  output  electrical  signal  at  each 
element  of  a  planar  array  of  point  sources.   The  output 
signal  is  described  in  terms  of  the  frequency  spectrum  of 
the  transmitted  electrical  signal,  the  transmit  and  receive 
planar  arrays,  and  the  random  ocean  medium  transfer  function. 
Vos  [Ref.  3]  used  these  results  to  develop  a  computer  program 
that  generates  time-samples  of  a  real  baseband  output 
electrical  signal  at  each  element  in  a  receive  planar  array 
as  a  function  of  a  variable  ocean  medium  sound  speed  profile, 
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planar  array  size,  array  far-field  beam  patterns  and  func- 
tional form  of  the  transmit  signal.   Ziomek  has  since  modi- 
fied this  program  to  generate  time-samples  of  the  complex 
envelope  of  real  bandpass  output  electrical  signals. 

The  research  documented  in  this  thesis  has  the  following 
objectives : 

develop  a  computer  simulation  of  a  correlator 
receiver  which  processes  the  output  electrical  signals 
generated  by  the  computer  program  developed  by  Ziomek 
and  Vos; 

apply  the  concept  of  model-based  signal  processing 
to  the  development  of  the  signal  processing  algorithm 
used  by  the  receiver; 

determine  the  effectiveness  of  the  approach  in  the 
detection  of  signals  from  a  planar  array  of  point 
source  elements. 

Since  the  effects  of  the  ocean  medium  on  the  signals 

processed  by  the  receiver  are  embodied  in  the  random  ocean 

medium  transfer  function,  the  basic  question  to  be  addressed 

may  be  stated  as  follows.   Can  a  priori  knowledge  of  the 

ocean  medium,  based  on  physical  principles  of  acoustic  wave 

propagation,  be  used  to  improve  the  detection  of  signals  in 

a  receiver,  and  to  what  degree? 

Ziomek1 s  use  of  linear  systems  theory  to  develop  a 

transfer  function  model  of  the  ocean  medium  immediately 

suggests  the  use  of  a  compensating  filter  at  the  array 

output  to  remove  the  undesirable  time  delays  due  to  system 

geometry  and  wave  propagation  effects.   This  filter  would 

ideally  cophase  the  signals  at  each  element  in  the  planar 

array,  resulting  in  maximum  signal  output  when  the  signals 
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are  added  together.   This  filter  is  implemented  in  the 
frequency  domain  through  the  use  of  Discrete  Fourier 
Transforms  (DFT) ,  and  the  approach  is  exactly  analogous  to 
the  FFT  beamforming  procedure  discussed  by  Ziomek  [Ref.  2: 
pp.  153-176] . 

In  the  frequency  domain,  the  time  delays  due  to  system 
geometry  and  wave  propagation  effects  are  represented  as 
phase  shifts  which  may  be  eliminated  if  known  a  priori.   The 
concept  of  model-based  signal  processing  is  applied  here  to 
obtain  the  proper  compensating  phase  shift  for  the  known 
system  geometry  and  wave  propagation  conditions. 

Section  II  describes  the  theory  used  to  develop  the 
receiver  model.   The  system  context  within  which  the 
receiver  operates  is  described  and  related  to  previous 
investigations.   A  functional  description  of  the  receiver 
is  shown,  and  each  of  the  major  functional  blocks  is  explained 
in  some  detail.   Finally,  a  statistical  description  of  the 
receiver's  performance  is  developed. 

The  computer  implementation  of  the  receiver  structure  is 
described  in  Section  III.   The  logical  flow  of  the  computer 
program  is  discussed  and  related  to  the  receiver  descrip- 
tion.  The  use  of  multiple  trials  to  estimate  the  probability 
of  detection  is  explained.   Each  of  the  major  subprograms 
is  characterized  in  terms  of  function  and  implementation. 
Verification  of  the  computer  simulation  is  discussed  last. 
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Section  IV  presents  the  data  obtained  from  the  simula- 
tion when  a  rectangular-envelope,  continuous  wave  (CW) 
pulse  or  a  rectangular-envelope,  linear-frequency-modulated, 
(LFM)  carrier  is  transmitted.   Receiver  performance  is 
described  by  plotting  the  probability  of  detection  (Pd) , 
for  a  given  probability  of  false  alarm  (Pfa) ,  as  a  function 
of  the  input  signal-to-noise  ratio  (SNR)  at  each  element 
in  the  receive  array.   Plots  are  provided  for  different  values 
of  Pfa,  and  show  the  relative  improvement  in  receiver  per- 
formance as  various  medium  and  wave  propagation  effects 
are  compensated  for  by  the  model-based,  signal  processing 
algorithm.   In  each  plot,  the  receiver  performance  predicted 
by  theory  when  all  array  element  output  signals  are  precisely 
cophased,  and  the  array  element  .input  noise  is  zero-mean, 
uncorrelated  and  Gaussian  is  shown  as  a  dashed  line.   The 
dashed  line  is  plotted  from  data  obtained  from  a  closed 
form  expression  relating  Pd  to  Pfa  as  a  function  of  array 
element  input  SNR,  and  is  superimposed  on  the  output  plots 
to  provide  a  baseline  for  judging  the  validity  of  the  receiver 
simulation  output  data. 

Conclusions  and  recommendations  are  discussed  in 
Section  V. 
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II.   THEORY  OF  THE  RECEIVER  MODEL 

A  receiver  operates  in  the  context  of  a  total  communi- 
cation system  consisting  of  a  signal  source  (transmitter) , 
a  signal  propagation  medium  (channel)  and  a  signal  sink 
(receiver) .   It  is  the  model  of  the  communication  channel 
that  is  of  initial  interest  since  the  signal  processing 
algorithm  will  depend  in  large  part  on  the  physics  des- 
cribing the  propagation  of  the  signal  through  the  channel. 

A.   OVERVIEW  OF  THE  COMMUNICATION  SYSTEM 

Ziomek's  model  (Ref.  2]  of  the  ocean  medium  is  described 
in  general  as  a  time-variant,  space-variant,  random  filter 
(transfer  function)  in  which  the  index  of  refraction,  or 
equivalently ,  the  speed  of  sound,  is  a  function  of  depth, 
and  includes  both  a  deterministic  and  a  random  component. 
However,  in  describing  the  electrical  output  signals  from 
the  receive   aperture,  the  model  becomes  more  restrictive  in 
the  sense  that  the  channel  is  considered  to  be  time-invariant, 
but  still  space-variant.   Furthermore,  the  transmit  and 
receive  apertures  are  taken  to  be  rectangular,  planar  arrays 
whose  elements  consist  of  complex  weighted  point  sources. 
Complex  weighting  of  the  array  elements  provides  the  means 
for  amplitude  shading  and  beam  steering  both  the  transmit 
and  receive  array  patterns.   The  complex  weights  are  ideal 
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for  removing  the  undesired  effects  of  the  channel  on  the 
output  electrical  signals  from  the  receive  array  elements, 
and  become  the  tool  for  applying  the  model-based,  signal 
processing  concept.   Figure  1  depicts  the  geometry  of  the 
transmit  and  receive  arrays  in  the  transmission  medium. 
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Figure  1.   System  Geometry 

Ziomek  goes  further  than  simply  considering  system 
geometry  in  generating  the  output  signals  at  the  receive 
array.   The  spatial  variance  of  the  transfer  function  per- 
mits modeling  of  propagation  effects  due  to  the  dependence 
of  the  index  of  refraction  on  depth.   Figure  2  shows  the 
ray-bending  effects  of  propagation  through  such  an 
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Figure  2.   Ray  Path  Bending  Due  to  Inhomogeneous  Medium 

inhomogeneous  medium.   The  phase  shifts  due  to  system  geometry 
and  propagation  effects  may  be  considered  together  or 
separately  in  the  generation  of  the  output  electrical  signal 
data  in  the  computer  simulation  due  to  Vos  [Ref .  3] .   It 
is  this  output  electrical  signal  data  which  is  used  as  the 
input  signal  to  the  receiver. 

B.   FUNCTIONAL  DESCRIPTION  OF  THE  RECEIVER 

The  output  electrical  signal  data  processed  by  the 
receiver  are  time  samples  of  the  complex  envelope  of  the 
bandpass  electrical  signal  at  the  output  of  each  element  in 
the  receive  array.   This  implies  that  the  bandpass  acoustic 
signal  incident  on  the  receive  array  has  been  converted  to 
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an  output  electrical  signal  by  each  transducer  element  in 
the  array.   Each  real,  bandpass  electrical  signal  then 
passes  through  a  quadrature  demodulator  to  become  an 
equivalent  baseband,  complex  envelope  signal,  which  is 
time-sampled  and  converted  from  analog  to  digital  form.   The 
complex  envelope  is  represented  by  the  I-channel  (in-phase) 
and  Q-channel  (quadrature-phase)  components  generated  by  the 
quadrature  demodulator.   Time-sampling  is  done  in  a  manner 
that  satisfies  the  Nyquist  criterion  for  the  baseband  infor- 
mation contained  in  the  I  and  Q  channels.   Thus,  many  of 
the  components  associated  with  a  receiving  system  are  already 
contained  within  the  simulation  that  generates  the  planar 
array  output  signal  data.   The  receiver  simulation  assumes 
these  components  exist,  and  essentially  provides  signal 
processing  of  the  complex  envelope  of  the  output  bandpass 
electrical  signals.   The  major  functional  blocks  of  the 
receiver  include: 

-  an  array  signal  processor, 

-  a  correlator  implementation  of  a  matched  filter 
receiver, 

-  a  magnitude-square  operation,  and 

-  a  threshold  decision  operation. 

Figure  3  shows  the  functional  block  diagram  for  the  receiver 
model.   The  magnitude-square  and  the  threshold  decision 
operations  are  considered  to  be  part  of  the  correlator 
(matched  filter)  detector  function  block  shown  in  Figure  3. 
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Figure  3.   Receiver  Block  Diagram 

The  receive  array  is  assumed  to  be  a  rectangular  planar 
array  of  M  xN  elements  where  both  M  and  N  are  odd  numbers. 
Each  element  is  assumed  to  be  an  omnidirectional  point 
source.   The  geometry  of  the  planar  array  and  associated 
mathematical  notation  is  shown  in  Figure  4. 
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Figure  4.   Planar  Array  Geometry 
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The  function  of  each  quadrature  demodulator  is  to 
convert  the  amplitude  and  angle-modulated,  bandpass, 
electrical  signal  at  each  transducer  output  into  its 
baseband,  complex  envelope.   Thus,  the  low-pass  complex 
envelope  may  be  sampled  at  a  much  lower  rate.   Note  that  the 
use  of  a  quadrature  demodulator  does  imply  that  the  carrier 
frequency  of  the  transmitted  signal  is  known.   It  should  be 
emphasized  that  the  output  electrical  signal  at  each  element 
in  the  receive  array  is  passed  through  its  own  quadrature 
demodulator  before  array  processing  begins.   The  quadrature 
demodulator  is  shown  schematically  in  Figure  5,  where  the 
complex  envelope  of  y(t) ,  denoted  y(t) ,  is  given  by, 


y(t)       =      yc(t)     +    j  ys(t) 


(2.1) 


ARRAY 

TRANSDUCER 

ELEMENT 


2Cos(2nfct) 
2SJn  (27Tfct) 


I 


w  +w 

LOW  PASS  FILTER 


yc(t) 


Q 
-O 


y*(t) 


Figure  5.   Array  Element  Quadrature  Demodulator 
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where : 


yc(t) 


ys(t) 


is  the  I-channel  (in-phase)  component 
of  y ( t) ,  and 

is  the  Q-channel  (quadrature-phase) 
component  of  y(t). 


The  array  signal  processor,  shown  schematically  in 
Figure  6,  is  a  FFT  beamformer.   The  function  of  the  array 
processor  is  to  maximize  the  total  output  signal  when  the 
signals  from  each  element  in  the  array  are  summed. 


y(lTs,mdx,ndy)        Yr(qf0,mdx)ndy) 


(M> 


DFT 

w.  r.  t. 

TIME 

c       -a     eJUmb    eJVn 


YT(ITS) 


COMPLEX 

WEIGHT 

GENERATOR 

Figure  6.   Array  Processor  Block  Diagram 

The  array  processor  is  essentially  a  filter  that  compen- 
sates for  system  geometry  and  wave  propagation  effects  on 
the  signal  transmitted  through  the  channel.   The  array 
processor  is  implemented  in  the  frequency  domain  where 
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filtering  can  be  obtained  by  multiplication  of  each  spectral 
component  of  the  complex  envelope  by  a  complex  weighting 
coefficient.   The  complex  coefficients  are  computed  from 
a  knowledge  of  the  channel  model,  thus  the  application  of 
model-based  signal  processing.   After  applying  the  complex 
weights,  the  Inverse  DFT  is  computed  to  recover  cophased, 
time  domain,  complex  envelope  signals  from  each  array  ele- 
ment.  By  summing  these  cophased  signals  over  all  MxM  array 
elements,  a  constructive  interference  effect  is  achieved, 
and  the  total  output  signal  is  maximized. 

The  matched  filter  portion  of  the  receiver  is  imple- 
mented by  correlating  the  total  time-sampled  output  signal 

and  noise  from  the  array,  r(£T  ).   Where, 

s 


r(£T  )   =   y^UT  )  +  n_UT  )  (2.2) 

s       is     is 


with  a  time  and  frequency  shifted  replica  of  the  complex 
envelope  of  the  transmitted  waveform,  x(t).   Since  the 
phase  of  the  received  signal  is,  in  general,  unknown,  the 
magnitude-square  of  the  correlator  output  is  taken  as  the 
input  to  the  threshold  detector.   This  input  is  compared  to 
a  preset  threshold  level  y    to  determine  the  presence  of  a 
signal.   The  preset  threshold  y    is  computed  from  a  Neyman- 
Pearson  criterion.   The  schematic  of  the  correlator/matched 
filter  detector  is  shown  in  Figure  7,  and  has  been  shown  by 
Van  Trees  [Ref.  4:pp.  244-247]  to  be  the  optimum  receiver 
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Figure  7.   Correlator/Matched  Filter  Detector 

for  the  detection  of  a  bandpass  signal,  with  random  ampli- 
tude and  phase,  in  the  presence  of  white,  Gaussian  noise. 

C.   STATISTICAL  DESCRIPTION  OF  THE  RECEIVER 

The  output  electrical  signal  data  embodies  several 
assumptions  and  restrictions  that  must  be  restated  before 
further  development  of  the  receiver  model.   The  signal  data 
is  assumed  to  be  generated  by  a  real,  bandpass,  acoustic 
field,  y.-(t,r) ,  incident  upon  an  array  of  electroacoustic 
transducers.   The  acoustic  wave,  y  (t,r),  is  propagating 

in  the  +n   direction  with  velocity  c,  or 

o  J 


yM(t,r)   =   x(t  +  [£-no]/c) 


(2.3) 


where ; 
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/\         /\ 


r   =   xx  +  yy  +  zz  (2.4) 


n    =   ux+vy+wz  (2.5) 

o      o     o     o 


u    =   sin  6   cos  \b  (2.6) 

o  o      o 


v    =   sin  6  sin  \b  (2.7) 

o  o     o 


and 


w    =   cos  B  (2.8) 

o  o 


Note  that  x(t)  may  be  an  arbitrary  function  of  time. 
1 .   Array  Element  Output  Signal  Description 

The  acoustic  field,  incident  upon  an  element  in  the 
array,  is  converted  to  an  output  electrical  signal  by  the 
transducer,  and  is  transformed  into  a  baseband  complex 
envelope  at  the  output  of  the  quadrature  demodulator.   After 
time  sampling  at  a  rate  f   =  1/T   (in  samples  per  second) , 
spatial  sampling  over  the  receive  planar  aperture  in  incre- 
ments of  d   and  d  ,  and  assuming  linearity  in  the  transducer 
x       y  3  J 

operation,  the  complex  envelope  of  the  output  electrical 
signal  at  the  m   ,n    element  of  the  planar  array  lying  in 
the  XY-plane  may  be  written  as, 


yUTs,mdx,ndy)   =   xUTg  +  [uQmdx  +  vQnd  ]/c)  (2.9) 
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For  example,  if  the  real,  bandpass  output  electrical 
signal  y(t)  is  an  amplitude  and  angle  modulated  cosine  wave 
of  the  form 


y(t)   =   a(t)cos[27Tf  t  +  6  (t)  ]  (2.10) 


after  quadrature  demodulation,  the  complex  envelope  of  y(t) 
may  be  represented  as 


y(t)   =   yc(t)  +  jys(t)  (2.11) 


where : 


y  (t)   =   a(t)cos  0(t)  (2.12) 

■c 


and 


y  (t)   =   a(t)sin  6(t)  (2.13) 

-1  s 


or,  in  magnitude-phase  form 


t)   =   ^yl(t)    +  vlit] 


y"(t)  +  yc(t)  (2.14) 


and 


(t)   =   Tan  1[ys(t)/yc(t)]  (2.15) 
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It  is  also  possible  to  show  that  the  envelope  function, 
E(t),  of  the  real,  bandpass  signal  y(t)  is 

E(t)   =   |a(t) |  (2.16) 

The  main  advantage  of  working  with  the  complex 
envelope  form  of  the  real  bandpass  signal  is  that  the 
baseband  complex  envelope  may  be  sampled  at  a  rate  f   (in 
samples  per  second)  determined  by  the  bandwidth  of  the 
baseband  modulating  waveform,  independent  of  the  carrier 
frequency . 

It  is  further  assumed  that  the  complex  envelope  of 
the  transmitted  waveform  x(t)  can  be  represented  exactly, 
over  an  interval  of  T   seconds,  by  a  finite,  complex, 
Fourier  series  such  as, 


K 
x(t)   =    I       0     exp  [  +  j27rqf  t]  (2.17) 

q=-K   q 


where  the  complex  Fourier  coefficient  c   can  be  written  as, 

q 


c    =   a  exp  [+j6  ]  (2.18) 

q     q      q 


where 


f  :   is  the  fundamental  frequency  in  Hertz  of 
the  signal  x(t)  with  period  T   =  l/fQ 
seconds ,  and 


26 


K:   is  the  maximum  number ^of  harmonics  used  to 
represent  the  signal  x(t) . 


The  complex  Fourier  coefficients  c   can  then  be 

q 

determined  directly  from  the  DFT  with  respect  to  time  of 
x(£T  ) ,  or 


T  • 

0 

cM   =  T        J  x(£Tc)W^  (2.19) 

£=-L 


q  L    0__T  •  S     L 


where 


WT   =   exp[  +  j2-^/L]  (2.20) 

J_l 


L'   =   (L-D/2  (2.21) 


and 


L:  is  the  total  number  of  time  samples  taken 
during  the  time  interval  T  seconds  where 
L  is  a  non-negative,  odd  integer, 


T0 :   is  the  fundamental  period  or  data  record 
length  in  seconds,  and 


o 


T  :   is  the  sampling  period  in  seconds  (note 
fs  =  1/TS) • 

To  satisfy  the  Nyquist  sampling  theorem,  it  can  be 
shown  [Ref.  2:pp.  164-165]  that, 

L   >   2K  +  1  (2.22) 


27 


and, 


T    <   T  /L  (2.23) 

s   —   o/ 


2 .   Generation  of  Complex  Weight  Phase  Factors 

With  these  assumptions  it  can  be  shown  [Ref.  2:pp. 
165-166]  that  the  normalized  DFT  with  respect  to  time  of  the 
complex  envelope  of  the  output  electrical  signal  y(£T  ,md  , 
nd  )  can  be  expressed  as, 


Y^  (qf  ,md  ,  nd  )   =   c  exp[  +  je  ]exp[xj2Tifu  m .,  /c] 
SN    o    x    y        q        m    r  tJ     o  dx/ 


where : 


and : 


x  exp  [  +  jc|)n]exp[^j27rfvond   /cj  (2.24) 


f      =      fc    +   qfQ    ;       -K, . . . ,q,. . . ,K  (2.25! 


f  :   is  the  carrier  frequency, 

c  :   is  the  constant  speed  of  sound  in  homogeneous 
ocean  medium  at  the  receive  array, 

u  :   is  the  direction  cosine  of  the  wave  propa- 
gation vector  along  the  x-axis,  equation 
(2.6) ,  and 

v  :   is  the  direction  cosine  of  the  wave  propaga- 
tion vector  along  the  y-axis.   See  equation 
(2.7) . 
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Phase  factors  0   and  <b      are  due  to  the  separable 

m      Tn  ^ 

complex  weight,  c   .   That  is, 

r        3      mn 


c  „   =   c  d       a  exp[  +  j6  ]b  exp[  +  jd)  ]  (2.26) 

mn       mn       m   ^   J  m   n   ^   J  Tn 


Notice  that  by  a  proper  choice  of  6   and  <p      in  equation 


(2.24)  where, 


)    =   ±j2TTfu  md  /c  (2.27) 

m       J     o   x 


and 


)    =   ±j27rfv  nd  /c  (2.28) 

n       J     o   y 


the  phase  shifts  due  to  system  geometry  may  be  completely 
cancelled  leaving  only  the  complex  Fourier  coefficients  c 
at  each  element  in  the  planar  array.   It  should  be  noted 
that  the  phase  correction  factors  are  functions  of  both  the 
array  element  position  (mdx  or  ndy) ,  and  the  frequency  f  of 
each  spectral  component  of  the  input  bandpass  signal  spec- 
trum.  That  is,  the  frequency  is  given  by  equation  (2.25). 

Once  the  phase  shifts  due  to  geometry  are  eliminated, 
taking  the  inverse  DFT  with  respect  to  frequency  of  the 
spectrum  of  the  output  electrical  signal  Y(qf  ,md  ,nd  ) 
will  yield  cophased  time  signals  at  each  element.   Summing 
all  the  cophased  signals  will  result  in  a  maximum  total 


29 


output  signal  y  (£T  )  from  the  array  processor.   That 
is, 


M '      N ' 
y  (£T  )   =     I  I  yUT  ,md  ,nd  )  (2.29) 

m=-M'  n=-N'  y 


where 


M1   =   (M-D/2  (2.30) 


and 


N'   =   (N-D/2  (2.31) 

This  approach  is  extended  to  compensate  for  the 
deterministic  signal  phase  shifts  caused  by  transmission 
through  an  inhomogeneous  ocean  medium.   Since  the  variation 
in  the  speed  of  sound  c(y)  is  assumed  to  be  a  function  of  y 
(depth)  alone,  the  calculation  involves  only  the  4>   phase 
factor  of  the  complex  weight.   Ziomek   has  shown  that  for  a 
sound-speed  profile  c(y)  with  a  constant  gradient  g,  a 
closed  form  expression  for  the  deterministic  component  of 
phase  shift  is, 


k   c 

*MD   =   +  2^T-{f  [nD(Y)  ~1J  +  AYn}  (2-32) 

o  ^ 


Extension  of  work  by  Ziomek  based  on  expression  for  6 
[Ref .  2:pp.  263-268] .  MD 
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where 


nD(y)   =   co/cD(y)  (2.33) 


CD(y)   =   cq  +  gAyn  (2.34) 


Ay     =    (y   -y  )   +  nd  (2.35) 

J  n       2  r  J  o  y 


kQ   =   2ufQ/c  (2.36) 


n  (y)  :   is  the  space-variant   (with  depth  only) 
index  of  refraction, 

cD(y) :   is  the  speed  of  sound  at  depth  y, 

c     :   is  the  speed  of  sound  at  the  transmit 
array, 

k     :   is  the  wave  propagation  constant  at 
the  transmit  array, 

g     :   is  the  gradient  (slope)  of  the  sound 
speed  profile, 

y     :   is  the  depth  of  the  center  element  of 
the  transmit  array,  and 

y     :   is  the  depth  of  the  center  element  of  the 
receive  array. 


The  negative  of  the  deterministic  medium  phase 
factor,  equation  (2.32),  is  simply  added  to  the  system 
geometry  phase  factor,  equation  (2.28) ,  to  obtain  the  total 
complex  weight  phase  factor  in  the  y  direction. 
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3 .   Array  Processor  Output  Signal  Statistics 

Development  of  a  receiver  model  also  requires  a 
specification  of  the  noise  environment  in  which  the  receiver 
operates.   Any  realistic  model  for  the  noise  environment  is 
extremely  complex.   However  in  order  to  form  tractable 
theoretical  results  that  can  be  reasonably  approximated  in 
a  computer  simulation,  zero-mean,  additive,  white,  Gaussian 
noise  (AWGN)  is  assumed  at  the  output  of  each  array  element. 
Remember  that  the  array  element  output  is  also  the  input  to 
the  quadrature  demodulator. 

The  AWGN  model  permits  derivation  of  a  closed  form 
expression  relating  Pd,  Pfa  and  array  input  Signal-to- 
Noise  Ratio  (SNR) .   This  type  of  noise  process  can  also  be 
reasonably  approximated  by  using  computer  generated  pseudo- 
random number  sequences  from  a  standard  Gaussian  random 
number  generator.   By  comparing  the  theoretically  predicted 
receiver  performance  with  the  results  of  the  simulation, 
verification  of  the  computer  implementation  of  the  receiver 
can  be  achieved.   Once  verified,  the  computer  simulation 
can  test  more  realistic  noise  models  with  some  confidence 
in  the  resulting  data. 

The  input  SNR  at  a  single  element  in  the  array  is 
defined  as, 


E{ |y(£T  ,md  ,nd  ) |  } 

SNRi   = - - ^ J"  (2.37) 

E{ |nUT  ,md  ,nd  ) |  } 
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where  E{  }  denotes  the  expected  value  of  the  quantity  within 
the  curly  braces. 

If  the  random  input  signal  y  ( £T  ,mdr,nd  )  is  ergodic, 

S     X     y 

then  the  mean-square  value  of  the  signal  can  be  found  by 
computing  the  time-average  instead  of  the  ensemble  average,, 
that  is, 


,2        ,  ~  ,2 

E{  y(£T  ,md  ,nd  )   }   =   <  y  (  £T  ,md  ,  nd    >  (2.38 

.  |Jsxy'         '  ■*        s   x   y  ' 


where  < ( • ) >  denotes  the  time-average  of  the  quantity  within 
the  parenthesis,  or 


T  /2 

l         ° 

<(•)>   =   sr-       /      C)  dt  (2.39) 

°  "To/2 

The  integral  in  equation  (2.39)  can  be  approximated 
for  computer  simulation  purposes  by, 


1    L' 
<(•)>   =  ^        I  (-)T  (2.40) 

o   =-L* 


where  the  infinitesimal  dt  is  approximated  by, 


dt   =   Ts  (2.41) 
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Because  of  the  assumption  regarding  the  use  of  a 
finite  Fourier  series  to  represent  the  original  transmitted 
waveform,  the  input  signal  power  at  a  single  element  in  the 
receive  array  can  be  computed  from  the  sum  of  the  magnitude- 
square  of  the  complex  Fourier  coefficients  [Ref.  5:pp.  44- 
45] ,  that  is, 


2         K       2 

<|yUT  ,md  ,nd  )  I  >  =     Y   |c  I  (2.42) 

1 2        s   x   y  '         '■  ,  '  q ' 
J  q=-K    ^ 


These  coefficients  are  easily  obtained  for  each  array 
element  by  computing  the  DFT  with  respect  to  time  of  the 
complex  envelope  of  the  output  electrical  signal  data  at 
each  element,  or 


L ' 

c    =   ^   7    y(£T  ,md  ,nd  )W?£  (2.43) 

q      L  isiL,  y         s    x'   y   L 


The  mean-square  value  of  the  zero  mean  noise  signal 

2 

at  the  array  element  input  is  equal  to  the  variance  o    of 


the  noise  input,  or 


E{|n(£T  ,md  ,nd  )|2}   =  q2  (2.44) 

1     s    x    y  '         mn 


Substituting  equations  (2.42)  and  (2.44)  into  equation  (2.37), 
the  expression  for  input  SNR  at  element  (m,n)  may  be  written 
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as, 


SNR.        7    |c  |2/o2  (2.45) 

1        L  _,  '  q '  '  ran 
q=-K    ^ 


or,  by  rearranging  the  variance,  equation  (2.45)  becomes 


a2         7    |c  |2/SNR.  (2.46 

mn        L  ,      '  q  '  '    l 
q=-K    ^ 


The  measured  input  signal  power  at  an  array  element 
and  a  desired  input  SNR  parameter  value  SNR.  can  be  used  to 
obtain  the  noise  power  (variance)  required  to  scale  the 
output  from  the  random  number  generator.   The  ability  to  set 
a  desired  input  SNR  value  is  necessary  in  order  to  test 
receiver  performance. 

Because  the  DFT  and  IDFT  are  linear  operations,  the 

noise  statistics  at  the  output  of  the  array  processor  are 

still  Gaussian,  and  uncorrelated  in  both  spatial  and 

temporal  coordinates.   The  total  noise  signal  nm ( £T  )  at 

J  T         s 

the   output   of    the    array   processor   may   be   written    as, 


M1  N1 

nTUTs)       =  I  I        nUTs,mdx,nd    )  (2.47) 


m=-M'  n=-N' 


2  . 

The  variance  of  the  total  noise  o_  is  equal  to  the  mean 

square  value  of  the  total  noise  signal  from  the  array 
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processor,  and  is  also  equal  to  the  sum  of  the  variances 
of  the  noise  signals  at  each  element  since  the  noise  process 
at  one  element  is  by  assumption  independent  of  the  noise 
process  at  any  other  element  in  the  array.   The  variance 
of  the  total  noise  signal  is  then 


E{|nUTs)|2)   =     J  I         omn  (2.48) 

1  m=-M'  n=-N' 


2 

If  the  variance  a  „  at  each  element  is  the  same  for  all 

mn 


elements  in  the  array,  or 


a2    =   a2  (2.49) 

mn 


2 

then,  the  variance  of  the  total  noise  a   can  be  written  in 

terms  of  the  variance  of  the  noise  at  each  element  in  the 
array  by  substituting  equation  (2.49)  into  equation  (2.48) 
to  obtain, 


a2   =   MNo2  (2.50) 


or,  in  terms  of  mean  square  values, 


E{  |nrT,(£Tc)  |2}   =   MNE{  |n(£To)  |2}  (2.51) 

1      o  S 
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i  2 
where  E{  n(£T  )   }  is  the  noise  power  at  an  array  element 

when  the  noise  power  is  considered  to  be  the  same  for  all 

array  elements  (m,n). 

Both  the  output  signal  and  output  noise  components 

at  all  elements  in  the  array  are  baseband  complex  envelope 

signals  with  I  and  Q  channel  components.   The  bandwidth  of 

the  noise  signal  is  set  by  the  bandwidth  of  the  low-pass 

filter  in  the  quadrature  demodulator  to  W  Hertz.   For  the 

purpose  of  the  simulation,  W  is  always  adjusted  to  include 

the  highest  harmonic  of  the  complex  envelope  signal,  that 

is, 


W   =   KfQ  (2.52) 


Knowing  the  bandwidth  W  and  using  the  AWGN  assumptions, 
the  noise  power  spectral  density  N   at  each  element  can  be 
related  to  the  variance  of  the  noise  process  at  each  element 
by, 


N    =   ai/(2Kf  )  (2.53) 

o       mn      o 


The  SNR  at  the  output  of  the  array  processor  SNR   is 
defined  as  the  ratio  of  the  total  signal  power  to  the  total 
noise  power,  or 


A   E{|y  (£T  ) |2} 

SNRa   =   ^ j-  (2.54) 

A      E(|n  UT  )  |Z} 
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where, 


2  M'     N1  2 

E{|yTUT  )|  }   =   E{|   I  I      yUT  ,md  ,nd  )  |  }      (2.55) 

1    S  m=-M'  n=-N'     s    x    y 


For  the  case  of  perfectly  cophased  signals  with  identical 
signal  power  at  each  element,  the  right-hand  side  of  equa- 
tion (2.55)  reduces  to 


E{ |yT(£Ts) |2}   =   (MN)2E{ |y(£Ts) |2}  (2.56) 


i  2 
where  E{[y(£T  ) |  }  is  equivalent  to  the  time-average  power 

i  2 
at  each  element  in  the  array,  and  E{|y  (£T  ) |  }  is  equivalent 

J.        O 

to  the  time-average  power  of  the  total  output  signal  from 
the  array  processor.   Using  the  uncorrelated  and  equal 
array  element  input  power  noise  assumptions,  the  expression 
for  SNR   can  be  rewritten  by  substituting  equations  (2.51) 
and  (2.56)  into  equation  (2.54),  or 


E{  |yUT  )  |  } 

SNR    =   MN  5—   =   MNSNR.  .  (2.57) 

A         E{|nUT  )  \     }  X 


Since  the  array  gain  (AG)  is  defined  as, 


AG   =   10  log   [SNRA/SNRd]   (dB)  (2.58) 
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the  array  gain  in  this  case  becomes, 


AG   =   10  log    [MN]   (dB)  (2.59 


Note  that  for  a  5  x 5  element  planar  array,  M  *N  =  25,  and 
the  array  gain  is  13.979  dB . 

4 .   Hypothesis  Testing  and  the  Neyman-Pearson  Criterion 
The  correlator/matched  filter  detector  portion  of 
the  receiver,  Figure  7,  is  modeled  as  a  binary  hypothesis 
testing  problem  using  the  Neyman-Pearson  decision  criterion. 
The  two  hypotheses,  H   and  H, ,  are  defined  as, 


I    V£V  +  nT 

U 


£T  )    :   H, 
s        1 

IT  ) 


T(£Ts)  :   Hq 


where,  in  terms  of  the  transmitted  waveform  x(t), 


M'    N1 
y  (IT  >   =   c    I  I         x  K«s-V),^,nd] 

m=-M'  n=-N'  2 


x  exp  [  +  j2TT(f)2,ilT  J  (2.61) 


A   s 


and  n  (£T  )  is  given  by  equation  (2.47).   If  one  assumes 

J.        O 

that  all  array  element  output  signals  are  cophased  and 
identical,  equation  (2.61)  reduces  to 
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y  (£T  )   =   cMNxUT  -t  )  exp  [  +  j  2  tt  4>  a  £T  ]  (2.62) 


where : 


and, 


'A 


c   =   a  exp  [+j6]  (2.63) 


t  :  is  the  actual  time  delay  in  seconds  at 
each  element  in  the  array  due  to  range 
separation  between  individual  elements 
of  the  transmit  and  receive  arrays, 

ta:    is  the  actual  time  delay  in  seconds  due 
to  range  separation  between  transmit 
and  receive  arrays  when  all  signals  are 
cophased, 


is  the  actual  doppler  shift  in  Hertz  due 
to  relative  motion  between  transmit  and 
receive  arrays, 


a:     is  the  amplitude  attenuation  factor  that 
is,  in  general,  a  random  variable, 

6:     is  the  generalized  phase  shift  of  the 

received  signal  with  respect  to  the  trans- 
mitted signal.   (In  general,  it  is  also  a 
random  variable  dependent  on  both  spatial 
and  temporal  coordinates.),  and 


Hq :    is  the  null  or  noise  only  (no  signal) 
hypothesis . 


Matched  filtering  is  obtained  by  correlating  the 
complex  envelope  of  the  total  array  output  r(£T  )  with  the 
complex  conjugate  of  the  processing  waveform  g(£T  ).   The 
functional  form  of  g ( £T  )  can  be  written  as, 


g(£T  )   =   x(£T  -T)exp[  +  j27Tcf)£T  J  (2.64) 

o  S  S 
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which  is  a  time  and  frequency  shifted  replica  of  the 
transmitted  signal  x(t)  where 


t :   is  the  estimate  of  time  delay  at  the  receiver, 
and 

$ :   is  the  estimate  of  the  doppler  shift  at  the 
receiver. 


For  the  purposes  of  this  study,  the  estimates  of 
time  delay  and  doppler  shift  are  assumed  to  be  precisely 
correct,  or 


x   =   t  -  t_   =   0.0  (2.65) 

A 


and 


$   =   <j>   -  <f>   =   0.0  (2.66) 


The  correlation,  or  inner  product  I   between  two 
functions  r(t)  and  g(t),  is  given  by 


I      =      <r(t),g(t)>   =     /   r(t)g*(t)dt         (2.67) 


which  is  approximated  in  the  simulation  using  the  trapezoidal 
rule  approximation  to  the  integral,  that  is, 


L'-l   r(£T  )g*UT  )  +  r[U+l)T  ]g*[U+l)T_) 
I      =    I = — T    (2.68) 

JL=-L'  S 
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The  magnitude-square  of  the  correlator  output  is 
taken  for  two  reasons.   First,  the  phase  of  the  carrier 
frequency  waveform  is,  in  general,  unknown.   As  the  phase 
of  the  received  carrier  varies  with  respect  to  the  phase  of 
the  quadrature  demodulator  local  oscillator  (LO)  signal,  the 
output  of  the  quadrature  demodulator  would  vary  from  a 
maximum  negative  to  a  maximum  positive  value  depending  on 
the  phase  difference  between  the  carrier  and  the  LO  which 
usually  is  taken  to  be  a  uniformly  distributed  random  varia- 
ble between  0  and  2tt  radians.   The  change  in  polarity  of 
the  quadrature  demodulator  output  would  propagate  through 
to  the  output  of  the  integrator  in  the  correlator/matched 
filter  detector,  Figure  7.   Taking  the  magnitude  of  the 
integrator  output  ensures  that  the  input  to  the  threshold 
comparator  will  always  be  non-negative  regardless  of  the 
phase  difference  between  the  carrier  and  LO  waveforms. 
Second,  when  the  array  element  noise  statistics  are  assumed 
to  be  Gaussian,  the  square  of  the  magnitude  of  the  integrator 
output  yields  an  input  to  the  threshold  comparator  which 
can  be  described  statistically  by  exponential  density 
functions  for  both  H,  and  H»  signal  hypotheses.   As  will 
be  shown  in  the  following  derivations,  the  exponential  den- 
sity functions  can  be  easily  integrated  to  obtain  a  closed 
form  expression  for  Pd  and  Pfa  in  terms  of  the  SNR  at  the 
input  to  the  threshold  comparator.   The  output  of  the 
magnitude-square  operation  is  the  sufficient  statistic 
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on  which  the  binary  decision  is  made.   That  is,  choose  H. 
if, 


I 


>       y 


(2.69 


or  choose  H,  if, 


(2.70) 


Assuming  that  the  total  noise  n  (t)  is  a  baseband 
Gaussian  process  of  bandwidth  W,  the  conditional  probability 
density  functions  (pdf 's)  of  the  magnitude-square  of  the 

correlator  output  with  and  without  a  signal  present  can  be 

2 
shown  to  be  exponential.    The  conditional  pdf ' s  are  given 

by 


P(U|2|H0)   =   l/(2a^)exp[-U|2/(2a^)] 


(2.71) 


and, 


P(  |£ 


V 


l/(2o^)exp[-|£|2/(2a2)] 


(2.72) 


For  zero  mean,  AWGN,  and  cophased,  equal  energy  (power) 
signals  at  each  array  element  output,  the  variances  of 


Derivations  for  the  pdf's,  expressions  for  Pd,  Pfa  and 
the  decision  threshold  y   were  provided  by  Prof.  L.  J.  Ziomek 
in  private  communication. 
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the  magnitude-square  of  the  correlator  output  can  be  shown 
to  be 


where 


o 


=   E{  \l 


N 


}/2 


MNN  EV2 


(2.73) 


and 


°1   =   E{|£s|2  +  KN|2)/2 


[  (MN)  E(a  }  X(t,4>)  I2  +  MNN  E~]/2 

o  X 


(2.74) 


I.    |  : 


is  the  magnitude-square  of  the  correla- 
tor output  when  the  input  to  the 
receive  array  is  noise  alone, 

is  the  magnitude  square  of  the  correla- 
tor output  when  the  input  to  receive 
array  consists  of  signal  alone, 


N 


o 


X 


E(a^}: 


is  the  power  spectral  density  level  of 
the  noise  signal  at  each  array  element 
output, 

is  the  energy  in  the  transmitted  signal 
which  for  simulation  purposes  is  defined 
to  be  equal  to  the  energy  E~  in  the  local 
processing  waveform  g(t), 

is  the  mean-square  value  of  the  amplitude 

attentuation  factor  (note:   E{a  }  =  a2 

only  when  deterministic  effects  are 
considered) ,  and 


X(x,<j>) 


is  the  magnitude-square  of  the  auto- 
ambiguity  function.   Note  that  |X((x,d>)  |2 
Et  |  X  (t  ,  4>)  I  z  ,  or  in  our  case  where  t  =  0 
and  cf)  =  0,  |X(0,0)|2  =  E?  since 


XN(0,0) |2  =  1. 


[Ref.  2:pp.  190-191] 
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To  calculate  the  threhold  y  for  a  desired  Pfa  and 


input  SNR,  the  integral  of  the  pdf  for  the  null,  or  H 


0 


hypothesis,  is  set  equal  to  the  desired  Pfa,  and  the  result- 
ing equation  is  solved  for  the  threshold.   The  magnitude- 
square  operation  results  in  particularly  simple  threshold 
equation  when  the  input  noise  is  assumed  to  be  white,  zero 
mean,  and  Gaussian.   Computing  the  Pfa  yields 


+00 


Pfa  /         P( \l I2    Hn)d|£ |2 


Y 


=      exp[-y/(2oQ) ] 


=      exp[-y/(MNN    E-~)  ] 

O      X 


(2.75) 


By  solving  equation  (2.75)  for  the  threshold  y, 


y   =   MNN  E~  in   [1/Pfa] 

U  2\. 


(2.76) 


Once  the  threshold  is  obtained,  the  correlator 
output  pdf  for  the  H-,  hypothesis  may  be  evaluated  to  obtain 
the  Pd. 


Pd   = 


/    p(K 


E±)d\l 


=      exp[-y/2a1] 


exp{-y/[  (MN)  2E{a}2  X(t,<J>)  |  2+MNN   E~}J 
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(2.77) 


Since  in  our  problem  t  =  0  and  §   =    0 , 


X(x,<t>)  I2   =   |X(0,0)  I2   =   E?  (2.78) 


x 


and  if  we  let 


E{a2}E~ 
SNR|^,2   =   MN( -)  (2.79) 


then, 


MNE{a2}E~ 

Pd   =   exp{-y/[MNN  E~  ( = +  1)]}         (2.80) 

r  O  X      N 

o 


and 


Pd   =   exp{-y/[MNN  E~(SNR,-|2  +  1)J}  (2.81) 

ox      \  Jt  \ 


Substitution  of  equation  (2.75)  into  equation  (2.81) 
gives  the  desired  closed  form  expression  relating  Pd,  Pfa 
and  the  SNR  of  the  magnitude-square  of  the  correlator  output, 
that  is, 


1/U+SNR,  t;  ,2] 
Pd   =   Pfa         '  '  (2.82) 


This  result  agrees  in  form  with  the  result  given  in  Van 
Trees  [Ref.  4:pp.  246-247]  for  a  similar  single  channel 
receiver  model.   Figure  8  graphically  depicts  the  relationship 
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2 


Figure  8.   Density  Functions  of  the  Magnitude- 
Sqaure  Correlator  Output 


between  the  conditional  pdf's,  the  decision  threshold  y,  the 
Pfa  and  the  Pd. 

The  correlator  output  SNR,  equation  (2.79),  can  be 
related  to  the  input  SNR  at  a  single  element  in  the  array 
through  the  array  gain,  and  a  factor  resulting  from  the 
slightly  different  definitions  of  input  SNR  and  correlator 
output  SNR.   In  terms  of  array  element  input  signal  energy 
E -  ,  the  input  SNR  given  in  equation  (2.45)  may  be  rewritten 
as, 


SNR. 

l 


1/T  )E~ 
° y 


mn 


E{a  }E~ 
x 

T  a  „ 
o  mn 


(2.83) 
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where , 


E~   =   E{a2}Ex  (2.84) 


is  the  average  received  energy  at  a  single  element  in  the 
array  due  to  the  transmitted  signal  x(t),  and  from  equation 

(2.53) 


7  2kN 

o  „   =   2kf  N    =   —-2.  (2.85) 

ran         o  o       T 

o 


The  magnitude-square  correlator  output  SNR  can  be  written 
as , 


SNR,^,2  =      _s  _-  (2.86) 


£l       Ed^l2) 


where, 


E{\1     \A}       =       (MN)2E{a2}E-  (2.87) 


and, 


E{|AJ2}   =   MNNE;  (2.88) 


Substitution  of  equations  (2.87)  and  (2.88)  into  equation 
(2.86)  gives 
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E{a2}E- 

SNR|~,2   =   MN  - (2.89) 

\t\  0 


By  substituting  equation  (2.85)  into  equation  (2.83)  and 
rearranging  terms,  equation  (2.83)  becomes 


E{a2}E- 

SNRi   =    2KN   X  (2-90 

o 


By  rearranging  terms  in  equation  (2.90)  and  substituting 
the  result  into  equation  (2.89),  the  desired  relationship 
expressing  the  magnitude-square  correlator  output  SNR  in 
terms  of  the  array  element  input  SNR  at  a  single  element 
in  the  array  is  obtained.   That  is, 


SNR, ^,2   =   MN2KSNR.  (2.91) 


Writing  both  sides  of  equation  (2.91)  in  dB  form  yields, 

SNR, ^ I  2  (dB)   =   10  1og10(MN)  +  10  1og(2K)  +  SNR.  (dB) 

(2  92 
AG  [ 

where  AG  is  the  array  gain  given  in  equation  (2.59). 

To  summarize,  Section  II. C  provides  the  equations 
needed  to  implement  a  computer  simulation  of  the  receiver 
structure  described  in  Section  II. B.   The  implementation  of 
the  computer  simulation  is  the  subject  of  Section  III. 
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Ill .   COMPUTER  SIMULATION  OF  THE  RECEIVER 

The  computer  program  RCVR  simulates  receiver  operation 
through  straightforward  application  of  the  equations  and  the 
concepts  developed  in  Section  II.   Written  in  FORTRAN,  the 
computer  program  RCVR  consists  of  a  top  level,  or  main 
program,  and  nine  subprograms.   The  computer  program  will 
be  explained  from  a  functional  viewpoint.   That  is,  the 
algorithms  used  to  implement  the  receiver  simulation  will  be 
related  to  the  theoretical  development  outlined  in  Section 
II,  but  translation  of  these  algorithms  into  FORTRAN  state- 
ments will  not  be  discussed.   The  main  program  will  be 
described  first.   The  description  of  the  main  program  will 
be  followed  by  a  detailed  discussion  of  each  subprogram.   In 
addition  to  explaining  the  computer  simulation,  the  methods 
used  to  validate  the  receiver  simulation  output  data  will 
be  presented. 

A.   TOP  LEVEL  PROGRAM  DESCRIPTION 

The  organization  and  logic  flow  of  the  top  level,  or 
main  program,  is  shown  in  Figures  9a  through  9e.   The 
functions  of  the  main  program  include: 

-  initializing  the  simulation  run-time  environment, 

-  invoking  subprograms  in  the  proper  sequence  to 
process  the  input  signal  data. 

-  providing  control  logic  and  noise  generation  algorithms 
needed  to  measure  receiver  performance,  and 
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Figure  9a.   Program  RCVR  Flowchart 
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Figure  9c.   Program  RCVR  Flowchart 
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Figure    9d.       Program   RCVR   Flowchart 
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-  selecting  the  form  of  the  simulation  output  data. 
Initialization  of  the  run-time  environment  involves  reading 
input  data  from  file  to  internal  storage,  generating  the 
local  processing  waveform,  measuring  the  time-average  input 
signal  power  at  each  array  element,  computing  the  baseband 
signal  bandwidth,  and  generating  the  complex  weights  used 
by  the  array  processor.   The  input  data  is  read  by  a  call 
to  the  subprogram  READY  which  returns  a  set  of  simulation 
parameters  in  COMMON  storage  and  the  complex  envelope 
electrical  signal  data  for  each  array  element  in  matrix 
form.   The  local  processing  waveform  is  obtained  by  a  call 
to  the  subprogram  SGNLGN.   The  time-average  signal  power 
at  each  element  (m,n)  in  the  planar  array  is  found  by  using  . 
the  fact  that  the  original  transmit  signal  was  synthesized 
from  a  finite  Fourier  series,  and  applying  equation  (2.43) 
to  obtain  the  complex  Fourier  coefficients,  c  .   Once  the 
complex  Fourier  coefficients  at  each  element  in  the  planar 
array  are  obtained,  the  time-average  power  is  computed  using 
equation  (2.42).   The  baseband  signal  bandwidth  W  is  com- 
puted using  equation  (2.52).   The  separable  complex  weights, 
cm  and  d  inequation  (2.26)  are  generated  in  two  stages. 
First,  the  real  amplitude  factors,  am  and  b   in  equation 
(2.26),  are  computed  by  a  call  to  subprogram  AMPWGT.   Second, 
the  real  phase  factors,  9   and  <p      in  equation  (2.26),  are 
obtained  by  a  call  to  subprogram  PHSWGT.   The  complex  weights 
are  then  computed  by  combining  the  amplitude  and  phase 
factors  as  shown  in  equation  (2.26). 
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The  first  signal  processing  step  is  to  generate  the 
total  array  output  signal  in  the  absence  of  additive  noise 
by  a  call  to  subprogram  ARYPRO.   Subprogram  ARYPRO  which 
uses  the  complex  weights  c   and  d   to  cophase  the  planar 
array  output  signal  returns  a  total  array  output  signals. 
The  total  array  output  signal  energy  is  computed  using 
equation  (2.40)  for  later  use  in  calculating  a  SNR  at  the 
output  of  the  array  processor.   The  ratio  of  array  output 
SNR  to  the  input  SNR  defines  the  array  gain,  and  provides  a 
check  on  the  validity  of  the  data  generated  by  the  array 
processor  algorithm. 

Receiver  performance  is  measured  by  computing  a  relative 
frequency  estimate  of  the  Pd  over  the  specified  range  of 
array  element  input  SNR  values  when  the  Pfa  is  a  known, 
time- invariant  parameter.   Using  the  Pfa  and  a  range  of  SNR 
values  specified  by  the  programmer  in  a  receiver  control 
data  block,  the  main  program  establishes  nine  input  SNR 
values  for  which  the  estimate  of  the  Pd  will  be  computed, 
and  determines  the  number  of  trials,  or  runs,  to  be  used  in 
computing  the  estimate. 

At  point  B  in  the  flowchart  of  Figure  9c,  the  main 
program  enters  a  loop  that  initializes  the  noise  source, 
performs  an  array  gain  calculation,  computes  the  detection 
threshold  and  then  enters  a  second,  inner  loop  where  the 
correlation  detection  is  done.   After  exiting  the  inner 
loop,  the  relative  frequency  estimate  of  the  Pd  is  computed. 
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The  outer  loop  increments  through  the  input  SNR  values 
determined  earlier,  and  exits  the  loop  when  the  specified 
maximum  input  SNR  value  is  reached. 

The  first  step  in  initializing  the  noise  generator  is 
to  determine  the  noise  variance  (power)  at  each  array  element. 
The  power  level,  or  variance,  of  the  noise  samples  for  a 
particular  element  in  the  array  is  a  function  of  the  time- 
average  signal  power  at  each  element  and  the  array  element 
input  SNR.   The  variance  can  be  found  using  equation  (2.46). 
By  controlling  the  power  level  of  the  noise  process  at 
each  array  element,  the  receiver  performance  can  be  measured 
over  a  range  of  array  element  input  SNR  values  regardless 
of  the  time-average  signal  power  level  at  a  particular 
array  element. 

The  variance  computed  in  this  manner  is  the  variance  of 
the  complex  envelope  baseband  noise  signal.   The  variance  of 
the  real  baseband  noise  signal  produced  by  the  noise  generator 
subprogram  AWGN  in  the  I  or  Q  channel  is  one-half  the  vari- 
ance of  the  complex  noise  signal  since  the  complex  envelope 
is  the  sum  of  the  two  independent  I  and  Q  channel  noise  sig- 
nals.  Therefore,  the  variance  of  the  complex  envelope  of 
the  noise  is  divided  in  half  prior  to  being  passed  to  the 
noise  generating  subprogram  as  a  scale  parameter. 

The  variance  of  the  total  complex  envelope  noise  signal 

2 
a   is  obtained  by  summing  all  the  noise  variances  of  the 

2 
array  element  complex  envelope  signals,  a   ,  as  indicated 
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by  equation  (2.48).   A  power  spectral  density  of  the  total 

2 
noise  signal  is  computed  using  equation  (2.53)  where  o   is 

2 
substituted  for  o 

mn 

The  SNR  at  the  array  processor  output  is  computed  using 
equation  (2.54)  where  the  time-average  power  of  the  total 
array  output  signal  is  substituted  for  the  mean  square 
ensemble  average,  and  the  mean  square  ensemble  average  of 
the  total  noise  is  taken  to  be  the  same  as  the  total  noise 
signal  variance.   The  array  gain  is  computed  using  equation 
(2.58),  and  is  held  in  internal  storage  for  later  output  in 
tabular  format. 

The  detection  threshold  y    is  computed  using  equation 

(2.76)  where  the  total  noise  power  spectral  density  is 

substituted  for  the  MNN   factor  and  the  energy  of  the  local 

o  ^J 

processing  waveform  E~  is  used  instead  of  the  transmit 

^  g 

signal  energy  E~ .   For  the  purpose  of  this  simulation  however, 
the  energy  in  the  local  processing  waveform  is  the  same  as 
the  energy  in  the  transmit  waveform.   That  is,  E~  and  E~ 
are  equal. 

The  second,  or  inner  loop,  of  the  main  program  begins 
at  point  D  in  Figure  9c.   Within  the  inner  loop,  the  complex 
envelope  noise  signals  are  generated  for  all  array  elements, 
the  total  noise  output  signal  is  generated  by  the  array 
processor  subprogram  ARYPRO,  the  total  signal  and  total 
noise  are  summed  and  correlated  with  the  local  processing 
waveform,  and  the  signal  detection  decision  is  made.   The 
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inner  loop  terminates  when  the  number  of  iterations  through 
the  loop  equals  the  number  of  trials  allowed  to  form  the 
estimate  of  Pd. 

The  time-sampled,  complex  envelope  array  element  noise 
data  is  obtained  by  repeated  calls  to  the  noise  generator 
subprogram  AWGN.   Each  call  to  AWGN  returns  a  properly 
scaled  pseudorandom  number  representing  one  sample  of  the 
noise  process  in  the  I  or  Q  channel  at  a  particular  array 
element  with  L  time  samples  taken  for  each  of  the  I  and  Q 
channels.   The  noise  signal  data  is  stored  in  the  same 
matrix  form  as  the  array  element  output  electrical  signal 
data. 

The  complex  envelope  noise  data  is  submitted  to  the 
array  processor  subprogram  for  processing  in  exactly  the 
same  manner  as  the  input  signal  data.   The  result  is  an 
array  total  noise  output  signal.   Processing  the  noise  alone 
in  this  manner  provides  some  gain  in  execution  speed  and 
provides  the  flexibility  to  estimate  the  probability  of 
false  alarm  directly,  if  desired.   The  DFT  and  IDFT  are 
linear  operations,  and  the  principle  of  superposition  holds. 
Therefore,  the  addition  of  the  total  noise  and  total  signal 
at  the  output  of  the  array  processor  is  equivalent  to  adding 
noise  to  the  signal  at  each  element  prior  to  passing  the 
data  to  the  array  processor  subprogram. 

Correlation  of  the  sum  of  the  total  signal  and  total 
noise  with  the  local  processing  waveform  is  accomplished  by 
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using  the  trapezoidal  rule  to  approximate  the  integral  of 
equation  (2.67).   The  magnitude-square  of  the  correlator 
output  is  obtained  by  taking  the  complex  product 


l\2      =      I    •  I*  (3.1) 


The  threshold  detection  portion  of  the  receiver  is 
implemented  by  comparing  the  output  of  the  magnitude-square 
operation  to  the  decision  threshold,  using  a  simple  IF-THEN- 
ELSE  binary  branch.   The  number  of  hits,  or  times  the  output 
exceeds  the  threshold,  are  counted  and  stored.   The  process 
of  generating  noise  samples,  and  making  a  hit  or  miss  deci- 
sion continues  through  a  large  number  of  trials.   Since  the 
correlation  is  done  using  the  H,  (signal  plus  noise)  hypothe- 
sis, Pd  can  be  directly  estimated  using  the  ratio 

Pd   =   HITS/TRIALS  (3.2) 

The  first  approach  taken  to  determine  the  minimum  number 
of  trials  required  to  estimate  Pd  was  based  on  equation 
(2.82).   After  rearranging  terms,  equation  (2.82)  can  be 
written  as 


Pfa   =   PdU+SNR]  (3.3) 


The  central  idea  was  to  compute  an  estimate  of  Pfa,  using 
equation  (3.3),  from  the  relative  frequency  estimate  of  Pd 
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in  equation  (3.2).   The  algorithm  would  then  terminate  when 
the  computed  estimate  of  Pfa  differed  from  the  Pfa  input 
parameter  by  some  arbitrary  small  amount.   However,  the  use 
of  equation  (3.3)  was  found  to  be  a  poor  test  for  establish- 
ing when  the  algorithm  should  terminate,  and  would  not 
terminate  the  algorithm  for  Pd  values  greater  than  about 
0.6. 

A  perturbation  sensitivity  analysis  of  equation  (3.3) 
can  be  performed  by  computing  the  total  differential  of 
equation  (3.3).   That  is, 


dpfa  ■  w  dPd  +  fit  dSNR 


dPfa   =   [l+SNRjPdSNRdPd  +  Pd [1+SNR:1  £n  (Pd)  dSNR      (3.4) 


Assuming  that  the  SNR  is  a  constant,  or  equivalently ,  the 
differential  dSNR  is  zero,  dPfa  becomes 


dPfa   =   [l+SNR]PdSNRdPd  (3.5) 


or 


APfa   =   [l+SNR]PdSNRAPd  (3.6) 


where 
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dPfa   *   APfa   =   |Pfa-Pfa|  (3.7) 


and 


dPd  z       APd      |Pd-Pd|  (3.8) 

Dividing  both  sides  of  equation  (3.6)  by  equation  (3.3) 
yields 


APfa.  =   [1+SNR]  APd 
Pfa  Pd 


where  APfa/Pfa  and  APd/Pd  represent  the  fractional  error 
between  the  actual  and  estimated  values  of  Pfa  and  Pd , 
respectively. 

Such  an  analysis  indicates  that  the  percent  error  in  the 
computed  estimate  of  Pfa  is  linearly  related  to  the  percent 
error  in  the  estimate  of  Pd.   The  constant  of  proportionality 
relating  the  error  in  the  computed  estimate  of  Pfa  to  the 
error  in  Pd  estimate  is  equal  to  (1  + SNR)  where  the  SNR  is 
taken  at  the  output  of  the  magnitude-square  operation.   Thus 
for  SNR  values  greater  than  approximately  1,  or  0  dB,  a  small 
error  in  the  relative  frequency  estimate  of  Pd  is  scaled  to 
a  larger  error  in  the  estimate  of  Pfa  found  by  using  equation 
(3.3).   Note  that  for  a  5  *5  element  planar  array  and  a 
bandwidth  of  5  times  the  fundamental  frequency  as  in  the  CW 
pulse  case,  the  array  element  input  SNR  can  be  found  by 
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equation  (2.92),  and  is  roughly  24  dB  less  than  the  SNR  at 
the  output  of  the  magnitude-square  operation.   That  is,  a 
0  dB  SNR  at  the  output  of  the  magnitude-square  operation 
corresponds  to  an  array  element  input  SNR  of  -24  dB .   The 
end  result  is  an  estimate  of  Pfa  that  diverges  wildly  from 
the  specified  Pfa  parameter  at  SNR  values  over  which  the 
receiver  operates.   In  addition,  equation  (3.3)  was  derived 
using  certain  assumptions  regarding  the  statistics  of  the 
input  noise  that  may  not  be  precisely  duplicated  by  the 
pseudorandom  noise  data  generated  by  the  computer  program. 
If  the  assumptions  regarding  the  use  of  zero  mean,  uncorre- 
cted, Gaussian  noise  are  not  satisfied  by  the  pseudorandom 
noise  source,  equation  (3.3)  does  not  hold,  and  the  algorithm 
would  not  be  suitable  for  other  input  noise  models.   For 
these  reasons,  the  use  of  equation  (3.3)  was  abandoned  in 
favor  of  an  empirically  determined  fixed  number  of  trials  to 
estimate  Pd. 

The  number  of  trials  needed  to  estimate  Pd  was  found  by 
assuming  that  the  absolute  minimum  number  of  trials  should 
be  the  reciprocal  of  the  Pfa  parameter.   That  is,  if  a  Pfa 
of  0.01  is  specified,  at  least  100  trials  must  be  taken  to 
allow  at  least  one  chance  in  one  hundred  of  a  false  alarm 
occurring.   Even  though  a  relative  frequency  estimate  of  Pd 
is  being  computed,  Pd  is  related  implicitly  to  Pfa  through 
the  decision  threshold,  and  the  value  of  the  Pfa  parameter 
should  be  taken  into  account  when  attempting  to  fix  the 
number  of  trials  needed  to  estimate  Pd. 
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The  number  of  trials  was  empirically  determined  by 
running  the  receiver  simulation  program  repeatedly  with 
increasing  multiples  of  the  minimum  number  of  trials,  and 
observing  the  effect  on  the  estimate  of  Pd.   Depending  upon 
the  value  of  the  Pfa  parameter,  it  was  found  that  four  to 
eight  times  the  minimum  number  of  trials  would  produce 
curves  that  did  not  change  appreciably  as  the  number  of 
trials  was  increased  further.   Therefore,  the  fixed  number 
of  trials  used  to  estimate  Pd  was  arbitrarily  set  at  10  x 1/Pfa 
for  Pfa  of  0.1,  and  to  5  x.l/Pfa  for  Pfa  of  0.01.   The  smaller 
multiplier  for  the  Pfa  of  0.01  became  necessary  due  to  limits 
on  computer  resources. 

1.   Subprogram  READY 

The  function  of  subprogram  READY  is  to  obtain  the 
simulation  parameters  and  the  complex  envelope  output 
electrical  signal  data  from  a  data  file.   The  data  file  is 
generated  by  the  ocean  communications  channel  simulation 
computer  program  developed  by  Vos  and  Ziomek  [Ref.  3].   That 
is,  READY  forms  the  interface  between  the  RCVR  simulation 
and  the  ocean  communications  channel  program.   A  flowchart 
of  READY  is  shown  in  Figure  10. 

The  simulation  parameters  are  read  in  first,  and 
are  stored  in  the  named  COMMON  blocks:   HEADER,  SIGNAL, 
ARRAY,  and  MEDIUM.   The  named  COMMON  blocks  provide  the 
mechanism  for  communicating  simulation  parameters  into  the 
various  subprograms.   The  HEADER  data  documents  the  type  of 
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Figure  10 .   Subprogram  READY  Flowchart 

communication  channel,  and  the  date  the  data  was  generated 
by  the  ocean  communication  channel  simulation.   The  SIGNAL 
data  includes:   the  fundamental  frequency  f0  and  the  period 
T0 ;  the  number  of  harmonics  K;  the  same  rate  f   and  the 
sample  period  T  ;  the  number  of  time  samples  L;  the  carrier 
frequency  f  ;  and  the  number  of  zeroes  padded  to  the  input 
signal  data,  if  any.   The  ARRAY  data  includes:   the  number 
array  elements  M  and  element  spacing  d   in  the  x-direction; 
the  number  of  elements  N  and  element  spacing  d   in  the  y- 
direction;  the  direction  cosines  u   and  v  ,  representing  the 
direction  of  the  direct  path  from  the  transmit  array  to  the 
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receive  array;  the  direction  cosines  u,  and  v,  ,  representing 
the  direction  the  transmit  array  beam  pattern  was  steered  in 
the  communication  channel  simulation;  the  depth  y   of  the 
center  element  of  the  transmit  array;  and  the  depth  y   of 
the  center  element  of  the  receive  array.   The  MEDIUM  data 
includes  the  speed  of  sound  c   at  the  center  element  of  the 
transmit  array;  the  sound-speed-profile  gradient  g;  and 
the  speed  of  sound  c    at  the  center  element  of  the  receive 
array. 

The  time-sampled,  complex  envelope,  output  electri- 
cal signal  data  is  read  in  next.   The  data  is  stored  in  the 
complex  matrix  variable  YCE  with  dimensions  L,  M  and 
N.   The  maximum  values  of  L,  M,  and  N  are  limited  to  33,  11 
and  11,  respectively. 

2.   Subprogram  SGNLGN 

The  function  of  subprogram  SGNLGN  is  to  generate  the 
time  samples  of  the  complex  envelope  and  to  compute  the 
energy  of  the  local  processing  waveform  given  by  equation 
(2.64).   A  flowchart  of  Subprogram  SGNLGN  is  shown  in  Figure 
11. 

The  local  processing  waveform  g(£T  )  is  synthesized 
from  a  finite,  complex  Fourier  series  with  provisions  for 
incorporating  the  estimates  of  time  delay  t  and  doppler 
shift  <J>  in  the  total  received  signal  y(ilT  ).   Accurate  esti- 
mates  of  T  and  <J)  are  necessary  to  achieve  maximum  output 
from  the  correlator/matched  filter  detector  portion  of  the 
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Figure  11.   Subprogram  SGNLGN  Flowchart 

receiver.   That  is,  the  maximum  receiver  sensitivity  is 
obtained  when  r  =  t   and  <J>  =  <J>   as  indicated  in  equations 
(2.65)  and  (2.66)  . 

For  the  purpose  of  this  study,  the  actual  doppler 
shift  is  always  set  to  zero  in  the  transmit  signal  syn- 
thesizer, and  the  actual  time  delay  due  to  the  range  or 
distance  between  array  is  set  by  the  system  geometry  under 
consideration  and'  the  reference  speed  of  sound  at  the 
transmit  array.   Therefore,  the  estimate  of  doppler  shift  is 
set  to  zero  in  subprogram  SGNLGN,  and  the  estimate  of  time 
delay  due  to  range  is  computed  from  the  line-of-sight  range 
between  the  center  element  of  the  transmit  array  and  the 
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center  element  of  the  receive  array,  and  the  speed  of  sound 
c   at  the  center  element  of  the  transmit  array,  or 


:x0-xr)2  +  (yQ-yr)2  +  (z0-zr)2/co     (3-10) 


where  (x  ,y  ,z  )  are  the  coordinates  of  the  center  element 
o   o   o 

of  the  transmit  array,  and  (x  ,y  , z  )  are  the  coordinates  of 
the  center  element  of  the  receive  array. 

The  estimates  of  x  and  $    are  easily  incorporated 
into  the  local  processing  waveform  of  equation  (2.64)  by 
applying  well-known  properties  of  Fourier  transforms  to 
equation  (2.17)  to  yield 


x(t-T)exp[j2Tr(J)t]   =    I     c   exp[j2irgf  (t-i) +j2irc|)t]      (3.11) 

q=-K  q  ° 


Note  that  the  right  hand  side  of  equation  (2.64)  is  just  the 
time-sampled  form  of  equation  (3.11). 

The  complex  Fourier  coefficients  c   of  equation 
(3.11)  are  identical  to  those  used  in  equation  (2.17)  to 
generate  the  transmit  signal  in  the  ocean  communication 
channel  simulation  computer  program.   Thus,  the  local 
processing  waveform  is  identical  to  the  transmit  signal  in 
functional  form  and  total  energy  content,  but  is  shifted  in 
time  and  frequency  by  the  estimates  of  range  delay  and 
doppler  shift,  respectively. 


69 


Subprogram  SGNLGN  also  computes  the  local  processing 
waveform  signal  energy  E~  for  later  use  in  setting  the 
decision  threshold  y.      To  compute  the  energy  the  magnitude- 
square  of  the  complex  Fourier  coefficients  are  summed  over 
all  harmonics  as  indicated  in  equation  (2.42)  where  c   is 
the  q-th  harmonic  of  the  local  processing  waveform.   This 
sum  is  equivalent  to  the  time-average  power  in  the  complex 
envelope  of  the  local  processing  waveform.   The  energy  can 
then  be  found  by  multiplying  the  time-average  power  by  the 
fundamental  period  T   of  the  local  processing  waveform. 

3.  Subprogram  AMPWGT 

The  function  of  subprogram  AMPWGT  is  to  provide  the 
real-valued,  amplitude  factors  am  and  b   of  the  separable 
complex  weights  c   and  d   in  equation  (2.26).   A  flowchart 
of  subprogram  AMPWGT  is  shown  in  Figure  12. 

To  generate  the  rectangular  amplitude  window,  a 
and  b   are  set  equal  to  1.0  for  all  elements  (m,n)  in  the 
receive  planar  array.   A  separate  subprogram  to  generate  the 
amplitude  weights  facilitates  generation  of  other  forms  of 
amplitude  windows  such  as  triangular,  Hamming,  Blackman, 
etc.   However,  only  the  rectangular  window  is  used  in  this 
study. 

4.  Subprogram  PHSWGT 

The  function  of  subprogram  PHSWGT  is  to  generate  the 
phase  factors  6   and  4>   of  the  separable  complex  weights 
c   and  d   in  equation  (2.26).   A  flowchart  of  subprogram 
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Figure  12.   Subprogram  AMPWGT  Flowchart 

PHSWGT  is  shown  in  Figures  13a  and  13b.   PSHWGT  can  be 
programmed  to  compute  phase  factors  that  compensate,  or 
remove  the  effects  of,  system  geometry  and  deterministic 
medium  wave  propagation  effects.   PHSWGT  can  also  introduce 
random  noise  in  the  phase  factors  for  study  purposes. 

The  phase  corrections  for  system  geometry  are  com- 
puted using  equations  (2.27)  and  (2.28),  where  the  direction 
cosines  u  and  v  are  selected  such  that  receive  array  beam 
pattern  is  aimed  at  the  transmit  array  along  the  direct  path 
from  the  receive  array  to  the  transmit  array.   That  is,  if 
n   in  equation  (2.5)  represents  the  direction  of  the  direct 
path  from  the  transmit  array  to  the  receive  array,  u  and  v 
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Figure  13a.   Subprogram  PHSWGT  Flowchart 
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Figure  13b.   Subprogram  PHSWGT  Flowchart 
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are  set  equal   to  -u   and  -v   so  that  the  receive  array 

o        o  J 

beam  pattern  points  in  the  -n   direction. 

The  deterministic  medium  phase  correction  factor 
is  found  by  using  equation  (2.32)  to  compute  the  phase  shift 
of  the  signal  due  to  propagation  through  an  inhomogeneous 
medium,  and  then  negating  the  result.   The  random  medium 
effect  can  also  be  computed  from  a  closed  form  expression, 
but  its  use  in  generating  phase  correction  factors  is  not 
the  subject  of  this  study,  and  will  not  be  discussed.   The 
total  phase  correction  factor  for  the  y-direction  is 
obtained  by  adding  the  system  geometry  and  deterministic 
correction  factors. 

5.   Subprogram  ARYPRO 

The  array  processor  subprogram  ARYPRO  uses  DFT  and 
IDFT  algorithms  to: 

-  generate  the  spectrum  of  the  input  electrical  signal 
data  at  each  element, 

-  correct  the  phase  of  the  spectral  components  to  co- 
phase  the  signals  at  all  array  elements,  and 

-  inverse  transform  the  signal  spectrum  at  each  array 
element  to  recover  the  cophased  signal  data. 

A  block  diagram  of  the  subprogram  ARYPRO  is  shown  in  Figures 

14a  and  14b. 

The  output  of  the  array  processor  is  a  time-sampled, 

complex  envelope  signal  representing  the  sum  over  all  array 

elements  of  the  signals  at  each  element.   The  defining 

equation  for  the  DFT  is  used  instead  of  a  fast  Fourier 

transform  algorithm  because  the  maximum  number  of  time 
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Figure    14a.       Subprogram   ARYPRO   Flowchart 
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Figure  14b.   Subprogram  ARYPRO  Flowchart 
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samples  permitted  (33)  is  insufficient  to  achieve  a  measurable 
improvement  in  speed  of  execution  [Ref.  6:pp.  151-152]. 
Additionally,  use  of  the  DFT  equation  relates  the  structure 
of  the  program  directly  to  previous  work  by  Ziomek  [Ref.  2] 
and  Vos  [Ref.  3] . 

6 .   Subprogram  AWGN 

The  function  of  subprogram  AWGN  is  to  generate  one 
sample  of  an  uncorrelated,  Gaussian  process  with  arbitrary 
mean  and  variance.   A  flowchart  of  AWGN  is  shown  in  Figure  15. 

AWGN  is  based  on  the  International  Mathematical 

Subroutine  Library  (IMSL)  FORTRAN  pseudorandom  number 

generator  routine  GGNQF.   GGNQF  is  a  function  subprogram 

that  returns  one  zero  mean,  unit  variance,  Gaussian,  or 

N(0,1),  pseudorandom  number  with  each  call  to  the  function 

subprogram.   The  zero  mean,  unit  variance  pseudorandom 

number  x   is  then  scaled  with  the  desired  mean  y   and  the 
n  x 

desired  standard  deviation  a      using  the  relation 


X  =   axXn  +  px  (3.12 


where  the  standard  deviation  is  computed  from  the  variance 
as 


a   =   \/a2  (3.13) 

x     V  x 


The  desired  mean  and  standard  deviation  are  passed 
to  AWGN  as  arguments  from  the  main  program.   The  mean  is 
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} 


Figure  15.   Subprogram  AWGN  Flowchart 

always  set  to  zero  for  the  purposes  of  this  study  and  the 
standard  deviation  is  found  from  equation  (3.13)  where  the 
variance  represents  the  "desired  power  level  in  either  the  I 
or  Q  channel  noise  signal. 

The  function  subprogram  GGNQF  internally  recomputes 
a  new  seed  value  for  subsequent  calls  from  a  seed  provided 
on  the  first  call  to  the  function.   Since  the  first  seed 
value  is  set  as  a  parameter  in  the  top  level  program,  and 
is  passed  as  an  argument  through  AWGN  to  GGNQF,  the  pseudo- 
random sequence  always  follows  the  same  pattern  each  time 
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the  receiver  simulation  is  run.   Using  the  same  pseudo- 
random sequence  for  each  simulation  run  is  essential  if 
comparisons  between  different  runs  are  to  be  made  from  the 
receiver  simulation  output  plots. 
7.   Subprogram  INTGRT 

The  function  of  subprogram  INTGRT  is  to  compute, 
using  equation  (2.68),  the  approximation  of  the  correlation 
integral  given  in  equation  (2.67).   A  flowchart  of  subprogram 
INTGRT  is  shown  in  Figure  16.   The  complex  product  of  the 
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Figure  16.   Subprogram  INTGRT  Flowchart 

total  signal  and  noise,  and  the  local  processing  waveform 
is  computed  in  the  top  level  program.   The  resulting 
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complex-valued  sequence  is  passed  to  INTGRT  which  computes 
the  correlator  output  I   by  using  a  trapezoidal  integration 
algorithm  to  approximate  the  integral  of  equation  (2.67). 
A  separate  subprogram  allows  easy  implementation  of  other 
numerical  integration  algorithms,  if  desired. 
8.   Subprogram  WRITBL 

The  function  of  subprogram  WRITBL  is  to  simply  pro- 
vide tabular  output  of  the  simulation  parameters  and  selected 
information  generated  during  the  execution  of  the  program. 
A  flowchart  of  subprogram  WRITBL  is  shown  in  Figure  17. 


WRITBL 


WRITE 

SIMULATION 

PARAMETERS 

TO  OUTPUT  FILE 


WRITE 

ESTIMATED 

PD  DATA 

TO  OUTPUT  FILE 


4 


RETURN 


) 


Figure  17.   Subprogram  WRITBL  Flowchart 

Besides  the  data  read  into  the  COMMON  storage  blocks 
by  subprogram  READY,  the  information  output  to  the  data  file 
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in  table  form  includes  the  array  element  input  SNR  at  which 
the  Pd  is  estimated,  the  estimate  of  Pd,  an  estimate  of  Pfa 
obtained  using  equation  (2.82),  the  SNR  at  the  output  of  the 
array  processor,  the  array  gain  and  the  SNR  at  the  output  of 
the  magnitude-square  operation.   The  tabular  output  is 
intended  primarily  to  provide  a  convenient  means  of  testing 
program  modifications.   It  is  not  intended  to  be  the  primary 
simulation  output,  and  will  not  be  discussed  further. 
9.   Subprogram  PDPLOT 

The  function  of  subprogram  PDPLOT  is  to  convert  the 
numeric  data  generated  by  the  receiver  simulation  into 
graphic  form.   A  flowchart  of  the  subprogram  PDPLOT  is  shown 
in  Figure  18.   The  plots  generated  from  the  numeric  data 
are  obtained  through  the  use  of  standard  DISSPLA  graphics 
library  subroutines. 

The  primary  output  of  the  receiver  simulation  is  a 
plot  of  the  estimate  of  Pd  versus  the  array  element  input 
SNR  with  a  fixed  value  of  Pfa  as  a  parameter.   PDPLOT  also 
computes  a  theoretical  value  of  Pd  using  the  equation  (2.82) 
where  the  SNR  at  the  output  of  the  magnitude-square  operation 
is  related  to  the  array  element  input  SNR  through  equation 
(2.91)  . 

The  curve  representing  the  estimated  Pd  is  plotted 
through  the  estimated  Pd  data  using  the  DISSPLA  least  squares, 
cubic  spline,  curve  fitting  plot  routine  SMOOTH.   The  value 
of  Pd  obtained  from  equation  (2.82)  is  plotted  using  the 
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Figure  18.   Subprogram  PDPLOT  Flowchart 
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DISSPLA  cubic  spline,  interpolating  polynomial  routine 
SPLINE,  and  appears  as  a  dashed  line  on  all  plots  generated 
by  the  receiver  simulation.   That  is,  the  curve  for  the 
calculated  value  of  Pd  is  an  interpolating  polynomial  that 
passes  through  the  calculated  Pd  data  points,  while  the 
curve  for  the  estimated  value  of  Pd  is  allowed  to  pass 
within  some  arbitrary,  small  offset  of  the  estimated  Pd 
data  point.   For  the  purpose  of  this  study,  the  maximum 
offset  allowed  was  10  percent  of  the  data  point  value. 
This  value  of  offset  provided  a  reasonable  compromise  be- 
tween obtaining  a  relatively  smooth  fit  of  the  plotted  curve 
without  diverging  excessively  from  the  data  points.   Both 
the  estimated  and  calculated  Pd  data  points  are  plotted  at 
the  same  array  element  input  SNR. 

B.   MODEL  VERIFICATION 

Verification  of  the  receiver  model  involved  two  tasks, 
characterization  of  the  pseudorandom  number  generator  noise 
source,  and  a  comparison  of  the  data  generated  by  the  re- 
ceiver simulation  program  with  the  results  predicted  by 
the  theoretical  development  discussed  in  Section  II.   Why 
simulate  a  test  case  for  which  a  theoretical  model  exists? 
The  main  reason  is  to  gain  some  confidence  in  the  results 
generated  by  the  simulation  when  the  inputs  are  such  that  a 
theoretical  model  does  not  exist  or  is  mathematically 
intractable. 


83 


1 .   Characterization  of  the  Noise  Source 

The  capability  to  generate  noise  samples  consistent 
with  the  noise  description  used  in  developing  the  receiver 
equations  is  central  to  verification  and  usefulness  of  the 
computer  simulation.   Thus,  the  first  step  is  to  ascertain 
the  statistical  properties  of  the  noise  source  used  in  the 
computer  program.   Recall  that  the  assumption  required  in 
developing  the  relations  between  Pd,  Pfa,  decision  threshold 
and  SNR  of  the  magnitude-square  of  the  correlator  output 
involved  the  use  of  zero  mean,  Gaussian  noise  which  is  uncorre- 
cted in  both  spatial  and  temporal  coordinates.   The  noise 
source  was  tested  by  taking  a  sequence  of  noise  samples  for 
both  I  and  Q  channels  in  manner  identical  to  that  used  in 
the  receiver  simulation.   That  is,  the  subprogram  AWGN 
was  embedded  in  a  test  program  to  ascertain  the  statistical 
properties  of  the  complex  envelope  of  the  noise  signal.   The 
following  tests  were  applied  to  the  I  and  Q  channel  sample 
sequences  to  verify  agreement  with  the  noise  process 
assumptions : 

-  sample  mean 

-  sample  variance 

-  histogram 

-  autocovariance,  and 

-  estimated  power  spectral  density. 

The  test  algorithms  were  obtained  from  standard  IMSL 
procedures . 
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Figures  19  and  20  show  the  results  of  the  sample 
mean,  variance  and  histogram  calculations  for  the  noise 
processes  in  the  I  and  Q  channels,  respectively.   The  sample 
mean  and  variance  calculations  approximate  the  zero  mean, 
unit  variance  assumptions  quite  well.   The  histogram 
calculations  indicate  the  distribution  of  the  samples  is 
roughly  Gaussian.   The  smooth  curve  plotted  through  the 
histogram  represents  the  exact  Gaussian  distribution  for  the 
size  of  the  sample  window  used.   The  greatest  deviation  from 
Gaussian  appears  near  the  mean  value  which  unfortunately  is 
where  most  of  the  sample  values  lie.   Thus,  some  difficulty 
in  achieving  a  perfect  correlation  between  simulated  per- 
formance and  theoretical  results  was  anticipated,  and  in 
fact  some  deviation  from  theory  at  large  values  of  Pd  did 
occur. 

The  space-time  correlation  properties  of  the  sample 
sequence  was  measured  by  computing  the  autocovariance  func- 
tion of  2000  samples  over  a  total  time  (or  space)  displace- 
ment (lag)  of  100  0  samples.   The  autocovariance  was  computed 
using  the  IMSL  routine  FTAUTO.   The  correlation  of  adjacent 
samples,  whether  one  or  two  or  1000  samples  apart,  was 
found  to  be  remarkably  small.   The  autocovariance  function 
of  the  I  and  Q  channel  noise  is  shown  in  Figures  21  and  22, 
respectively.   Because  of  the  manner  in  which  the  samples 
are  drawn  in  the  simulation,  the  distinction  of  a  sample 
being  assigned  to  a  particular  channel  at  a  particular 
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element  at  a  particular  time  is  immaterial.   The  auto- 
covariance  of  the  samples  indicates  that  noise  samples  can 
be  generated  which  are  independent  of  sample  time,  array 
element  location,  and  receiver  channel. 

The  estimated  PSD  was  computed  primarily  to  rein- 
force the  results  obtained  from  the  autocovariance  compu- 
tation.  The  IMSL  routine,  FTFPS ,  used  to  compute  the 
PSD  of  the  sample  sequence  implements  an  algorithm  similar 
to  that  due  to  Welch  [Ref.  7:pp.  553-554]  in  which  the  total 
sample  record  is  partitioned  into  contiguous  subrecords . 
Each  subrecord  is  amplitude  weighted  with  a  triangule  window 
function.   Then,  a  periodigram  of  each  amplitude  weighted 
sample  subrecord  is  computed.   The  resulting  subrecord 
periodigrams  are  averaged  over  all  subrecords  in  an  effort 
to  reduce  the  variance  of  the  estimate  of  the  PSD.   The 
estimate  of  PSD  for  the  I  and  Q  channel  noise  is  shown  in 
Figures  23  and  24,  respectively.   The  estimated  PSD  is 
approximately  flat  for  both  the  I  and  Q  channels  over  the 
sampled  frequency  range  of  0  to  tt  radians,  or  one  sample 
frequency  period.   This  tends  to  support  the  results  obtained 
from  the  autocovariance  computation  that  the  noise  samples 
are  indeed  uncorrelated. 

2 .   Verification  of  the  Output  Data 

The  primary  output  data  from  the  receiver  simulation 
is  a  plot  of  an  estimate  of  Pd  for  various  values  of  input 
SNR  at  a  given  Pfa.   As  shown  in  equation  (2.82)  ,  a  closed 
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form  relation  exists  between  the  Pd,  Pfa  and  correlator 
output  SNR.   Furthermore,  the  correlator  output  SNR  can  be 
related  to  the  array  input  SNR  when  the  signals  at  all  array 
elements  are  precisely  cophased  by  equation  (2.91). 
Therefore,  a  plot  of  the  theoretical  performance  computed 
using  equations  (2.82)  and  (2.91)  can  be  superimposed  on 
the  plot  of  the  receiver  simulated  performance.   The  curve 
representing  theoretical  performance  provides  verification 
of  the  simulated  performance  when  input  signal  parameters 
and  phase  weighting  allow  precise  signal  cophasing  to 
occur.   The  computed  performance  curve  also  provides  a 
baseline  to  evaluate  simulated  performance  when  the  input 
noise  or  receiver  operating  conditions  differ  from  the 
assumptions  used  in  formulating  the  receiver  model.   In 
all  the  plots,  the  theoretical  performance  is  shown  as  a 
dashed  curve. 


93 


IV.   PRESENTATION  OF  RESULTS 

This  section  presents  the  results  of  three  receiver 
simulation  case  studies  and  analyzes  the  simulation  output 
in  each  case.   The  three  cases  considered  are: 

-  Case  TEST,  in  which  the  output  electrical  signal  data 
at  each  array  element  is  produced  by  a  test  signal 
generator  computer  program,  and  is  free  of  communi- 
cation channel  effects  such  as  attenuation,  phase 
shifts  due  to  wave  front  refraction,  time  delay  due  to 
range  separation,  and  the  transmit  array  beam  pattern. 

-  Case  HMG1,  in  which  the  output  electrical  signal  data 
at  each  array  element  is  produced  by  the  ocean 
communication  channel  simulation  computer  program, 
and  includes  the  effects  of  path  attenuation,  phase 
shifts  due  to  system  geometry,  time  delay  due  to  range 
separation,  and  the  transmit  array  beam  pattern. 

-  Case  INHMG1,  in  which  the  output  electrical  signal 
data  at  each  array  element  is  produced  by  the  ocean 
communication  channel  simulation  computer  program, 
and  includes  not  only  the  effects  given  in  case  HMG1, 
but  also  contains  phase  shifts  due  to  refraction  of 
the  wavefront  along  the  propagation  path. 

The  results  of  these  case  studies  will  be  analyzed  by 
providing  a  brief  description  of  the  transmit  signal 
waveforms  used  to  generate  the  array  element  output  elec- 
trical signal  data,  by  listing  the  system  parameters  that 
distinguish  the  test  cases,  and  finally,  by  interpreting 
the  plots  of  Pd  versus  array  element  input  SNR  for  each 
case  study. 

A.   TRANSMIT  WAVEFORMS 

The  analysis  and  theoretical  development  of  the  receiver 
simulation  in  Section  II  makes  no  assumption  regarding  the 
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functional  form  of  the  transmit  signal.   The  receiver 
should  perform  equally  well  regardless  of  the  type  of  trans- 
mit signal  provided  that  the  correlator/matched  filter  uses 
a  replica  of  the  transmit  signal  as  the  local  processing 
waveform.   To  verify  this  hypothesis,  and  thereby  test  the 
integrity  of  the  receiver  simulation  and  the  validity  of  the 
assumptions  used  in  generating  the  model,  more  than  one  type 
of  transmit  signal  was  used  to  test  receiver  simulation 
performance.   The  use'  of  multiple  transmit  signal  waveforms 
also  provides  a  broader  base  from  which  to  draw  conclusions 
regarding  the  concept  of  model-based  signal  processing.   For 
the  purpose  of  this  research,  two  types  of  transmit  wave- 
forms were  used  in  producing  the  Pd  versus  array  element 
input  SNR  plots  in  each  of  the  three  case  studies. 

A  rectangular  envelope  CW  pulse  and  a  rectangular 
envelope  LFM  pulse  were  selected  as  transmit  signals.   These 
two  particular  forms  were  chosen  for  several  reasons.   First, 
the  CW  pulse  and  LFM  pulse  waveforms  were  considered  to  be 
typical  transmit  signals  used  in  SONAR  and  acoustic  signal- 
ing systems.   Second,  the  CW  pulse  provides  a  simple  case 
of  amplitude  modulation  while  the  LFM  pulse  provides  an 
example  of  an  amplitude  and  angle  modulated  waveform  with 
considerably  different  spectral  characteristics.   Finally, 
either  signal  can  be  readily  synthesized  from  a  finite, 
complex,  Fourier  series  whose  coefficients  c  may  be  calculated 
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from  a  closed  form  expression  for  the  frequency  spectrum 
of  the  time  domain  pulse  characteristics. 

The  rectangular  envelope  CW  pulse  is  discussed  first. 
The  pulse  repetition  frequency  is  the  same  as  the  funda- 
mental frequency  f   of  the  finite  frequency  spectrum  from 
which  the  pulse  is  synthesized.   The  pulse  duty  cycle  was 
arbitrarily  taken  to  be  0.5  yielding  a  pulse  width  of  T  /2 
seconds  where  T   =  1/f  .   The  complex,  Fourier  series 
coefficients  used  to  synthesize  the  complex  envelope  of  the 
CW  pulse  are  obtained  from  a  closed  form  expression  for  the 
complex-valued  continuous  spectrum.   To  obtain  the  Fourier 
coefficients,  the  closed  form  expression  for  the  continuous 
spectrum  is  evaluated  at  the  discrete  frequencies  qf  ,  and 
the  resulting  complex  value  is  divided  by  the  fundamental 
pulse  period  T  .   The  index  q  takes  on  integer  values 
between  -K  to  K  where  K  is  the  maximum  number  of  harmonics 
in  the  finite  Fourier  series  used  to  synthesize  the  CW 
pulse.   The  continuous  spectrum  of  the  CW  pulse  is  the 
familiar  sin(x)/x  form  obtained  by  taking  the  Fourier  trans- 
form of  the  rectangular  pulse  shape  where  the  pulse  width 
is  one-half  the  pulse  period.   The  following  specific 
transmit  signal  parameters  were  used  in  all  CW  pulse 
simulations : 

-  Amplitude,  A   =   40.0 

-  Duty  Cycle,  D    =   0.5 

-  Fundamental  Frequency,  f    =   200  Hz 
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-  Number  of  harmonics,  KMAX  =  5 

-  Harmonic  values.    c  =  20.00000  exp[j0°] 

c_1  =  c    =  12.33240  exp[j0°] 
c_„  =  c„   =  0.000000  exp[jO°J 
c_3  =  c3   =  4.244132  exp[jl80°] 
c_4  =  c.      =  0.000000  exp[j0°] 
c_5  =  c5   =  2.546479  exp[j0°] 
The  LFM  pulse,  or  linear  frequency  modulated  pulsed 
carrier  waveform,  is  discussed  next.   The  complex  Fourier 
coefficients  used  to  synthesize  the  LFM  pulse  are  found 
using  a  procedure  similar  to  that  used  for  the  CW  pulse 
except  the  closed  form  expression  for  the  complex-valued 
continuous  spectrum  of  the  LFM  pulse  was  found  by  using  the 
method  of  stationary  phase.   Officer  [Ref.  8:pp.  67-68] 
describes  the  method  of  stationary  phase  as  does  Papoulis 
[Ref.  9:pp.  267-273]  who  also  provides  a  complete  descrip- 
tion of  the  LFM  waveform.   The  following  transmit  signal 
parameters  are  used  in  all  LFM  pulse  simulations: 

-  Amplitude,  A     =  40.0 

-  Duty  Cycle,  D    =0.8 

-  Phase  Deviation  Constant        B     =  2356.2  radians/volt 

-  Fundamental  Frequency,  f     =  10  Hz 

-  Number  of  harmonics,  KMAX  =  3 

-  Harmonic  values,  c        =  14.60593  exp[j45°] 

c    =  c   =  14.60593  explj21°] 
c_2  =  c2  =  14.60593  exp[j309°] 
c    =  c   =  14.60593  exp[jl89°] 
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It  should  be  emphasized  that  these  complex,  Fourier 
series  coefficients  which  are  used  in  the  ocean  communica- 
tion channel  simulation  computer  program  to  synthesize  the 
complex  envelope  CW  and  LFM  pulse  transmit  waveforms  are 
also  used  in  subprogram  SGNLGN  to  produce  the  local  processing 
waveform  for  the  correlator/matched  filter  detector  section 
of  the  receiver  simulation  computer  program  RCVR. 

B.   CASE  TEST 

The  input  signal  data  for  case  TEST  is  obtained  from  a 
separate  computer  program  written  to  synthesize  array 
element  output  electrical  test  signals  from  a  finite  Fourier 
series  expansion  of  the  test  signal  waveform.   The  technique 
is  similar  to  that  used  in  generating  the  receiver  simula- 
tion local  processing  waveform  in  subprogram  SGNLGN.   How- 
ever, appropriate  phase  shifts  are  added  to  the  Fourier 
coefficients  to  vary  the  direction  of  arrival  of  the  plane 
wave  incident  on  the  receive  array,  and  to  incorporate  the 
time  delays  due  to  spatial  separation,  d   and  d  ,  of  the 
array  elements.   Furthermore,  the  time  delay  due  to  range 
and  the  doppler  shift  of  the  test  signals  are  both  set  to 
zero  to  prevent  signal  loss  in  the  correlator/matched  filter 
due  to  a  mismatch  between  the  total  array  output  signal  and 
the  local  processing  waveform.   In  this  way,  array  element 
output  electrical  signals  are  obtained  that  exclude  effects 
due  to  propagation  of  the  sound  wave  through  the  ocean 
medium,  or  spatial  modulation  of  the  signals  due  to  the 
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transmit  array  beam  pattern.   Since  the  medium  effects  are 
excluded  from  the  case  TEST  array  element  signals,  the 
carrier  frequency  f   is  set  to  zero,  and  the  array  element 
spacing  is  increased  to  keep  the  interelement  spacing  d 
and  d   equal  to  one-half  the  shortest  wavelength  in  the 
transmit  signal.   Note  that  this  produces  different  inter- 
element spacings  for  the  CW  and  LFM  pulse  situations. 

Signals  generated  in  this  fashion  allow  the  receiver 
output  to  be  analyzed  separately  from  the  communication 
channel.   Essentially,  it  is  equivalent  to  having  a  signal 
generator  that  provides  a  test  signal  to  measure  receiver 
performance.   The  results  obtained  from  case  TEST,  when 
combined  with  the  plot  of  theoretically  predicted  performance, 
provides  a  baseline  with  which  to  analyze  the  receiver 
output  when  the  ocean  communication  channel  simulation  data 
is  processed. 

The  system  parameters  defining  case  TEST  are  as  follows: 
-  Array  Parameters 

Number  of  array  elements,       M   =  N   =5 

Array  Element  Spacing,     CW    d   =  d   =  0.7375  meters 


LFM   d   =  d   =24.58  meters 
x     y 


-  Medium  Parameters 


Speed  of  Sound,  c   =  1475  meters/second 

Actual  time  delay  t   =  0.0  seconds 

Actual  doppler  shift  cf>A  =  0.0  Hertz 

Four  plots  representing  receiver  performance  were 
generated  for  case  TEST.   Phase  weighting  is  done  to 
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compensate  for  the  direction  of  arrival  of  an  acoustic  plane 
wave  incident  upon  the  array  elements.   The  direction  of 
arrival  can  be  set  to  any  arbitrary  direction  by  the  test 
signal  generator  computer  program.   However,  the  direction 
was  chosen  to  be  the  same  as  that  due  to  the  system  geometry 
used  in  the  case  studies  HMG1  and  INHMG1 .   Phase  weighting 
for  geometry  is  indicated  by  the  state  of  the  logical 
variable  STEER.   The  estimates  of  time  delay  due  to  range 
and  doppler  shift  were  both  set  to  zero  in  subprogram  SGNLGN 
to  maximize  the  output  of  the  correlator/matched  filter 
detector.   Two  of  the  plots,  Figures  25  and  26,  show  receiver 
performance  for  a  rectangular  envelope  CW  pulse.   Figure  25 
was  produced  using  a  Pfa  of  0.1,  and  Figure  26  shows  the 
effect  of  a  Pfa  of  0.01.   The  remaining  two  plots  generated 
for  case  TEST,  Figures  2  7  and  28,  show  receiver  performance 
for  a  rectangular  envelope  LFM  pulse  at  a  Pfa  of  0.1  and 
0.01,  respectively.   Since  the  input  electrical  signal  data 
does  not  incorporate  channel  effects,  phase  weighting  for 
deterministic  and  random  medium  effects  was  not  done  for 
case  TEST. 

Several  observations  regarding  model  performance  can 
be  made  by  analyzing  the  plots.   Generally  for  all  plots,  the 
measured  performance  of  the  simulation  (solid  curve)  follows 
the  predicted  theoretical  performance   (dashed  curve) ,  but 
the  agreement  is  not  exact.   That  is,  the  estimate  of  Pd 
follows  the  same  trend  as  the  predicted  Pd,  but  the 
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transition  from  a  low  Pd  to  a  Pd  of  1  is  more  sensitive  to 
changes  in  the  array  element  input  SNR.   A  comparison 
between  the  plots  shows  the  Pd  to  be  a  function  of  the  Pfa 
as  predicted  by  equation  (2.82),  and  for  small  values  of 
SNR,  the  estimated  Pd  approaches  the  Pfa  parameter  used  by 
the  simulation  in  generating  the  plots.   This  effect  is 
also  predicted  by  equation  (2.82)  since  at  SNR  values  much 
less  than  one,  the  input  signal  is  essentially  off,  and  only 
the  noise  is  producing  signal  detections.   A  comparison 
between  the  plots  for  the  CW  pulse  and  LFM  pulse  does  show 
some  minor  difference  in  receiver  performance,  but  the  same 
general  trends  seem  to  hold  regardless  of  the  functional 
form  of  the  transmit  signal. 

The  differences  in  the  measured  performance  between  the 
CW  pulse  and  the  LFM  pulse  simulations  can  be  attributed  in 
part  to  the  difference  in  the  maximum  number  of  harmonics 
used  to  synthesize  the  transmit  signals.   As  shown  by 
equation  (2.22),  the  number  of  time  samples  L  is  related  to 
the  maximum  number  of  harmonics  K  used  to  synthesize  the 
transmit  signal.   The  difference  in  the  number  of  time 
samples  causes  different  subsequences  of  noise  samples  to  be 
drawn  from  the  noise  generator  subprogram  AWGN.   The 
different  noise  sample  subsequences  will  cause  small  varia- 
tions in  the  relative  frequency  estimate  of  Pd.   Thus,  an 
exact  comparison  of  the  noise  performance  cannot  be  made 
across  different  types  of  transmit  signals  when  the  number 
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of  harmonics  used  to  create  the  transmit  signal  are  not  the 
same.   However,  the  real  significance  of  this  effect  is  that 
even  with  different  noise  sample  subsequences  in  the  I  and 
Q  receiver  channels  at  the  array  elements,  the  measured 
performance  conforms  closely  to  the  performance  predicted  by 
theory  under  the  very  ideal  assumptions  of  precise  signal 
cophasing,  and  zero-mean,  white,  Gaussian  input  noise. 

For  these  reasons,  it  can  be  inferred  that  the  basic 
receiver  model  is  valid,  but  some  difficulty  may  exist  in 
generating  precisely  white,  zero-mean,  Gaussian  noise  within 
the  simulation  program.   Indeed  as  indicated  in  Section 
III.B.l,  the  noise  source  used  in  the  receiver  simulation 
does  not  generate  precisely  zero-mean,  white,  Gaussian 
noise.   Furthermore,  the  measured  Pd  is  a  relative  frequency 
estimate  of  the  actual  receiver  performance,  and  may  also 
contribute  some  error  in  the  simulation  results. 

Considering  these  approximations  to  the  assumptions  used 
in  the  receiver  model  development  of  Section  II,  exact 
agreement   should  not  be  expected  even  in  the  case  where  the 
communication  channel  effects  are  excluded.   However  if 
only  relative  changes  or  effects  are  to  be  observed,  moderate 
disagreement  between  the  theory  and  simulated  receiver 
performance  can  be  negated  by  comparing  only  differences  in 
the  estimate  of  Pd  for  each  of  the  case  studies.   A  com- 
parison of  the  differences  in  the  measured  receiver  perfor- 
mance will  indicate  if,  and  to  what  degree,  the  effects  due 
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to  the  physics  of  acoustic  wave  propagation  can  be  compen- 
sated for  by  model-based  signal  processing.   The  results  of 
case  TEST  then  can  be  viewed  as  a  validation  of  the  receiver 
model  and  the  applied  programming  used  to  implement  the 
model. 

C.   CASE  HMG1 

The  input  signal  data  for  case  HMG1  and  case  INHMG1  is 
produced  by  the  ocean  communication  channel  simulation 
program  due  to  Vos  and  Ziomek  [Ref.  3].   For  case  HMG1 ,  the 
ocean  is  considered  to  be  a  homogeneous  medium  with  respect 
to  the  speed  of  sound  propagation  through  the  water.   That 
is,  the  speed  of  sound  is  identical  at  the  transmit  and 
receive  arrays,  and  at  all  points  in  between.   In  this  model 
of  the  ocean  medium,  refraction  or  ray  bending  of  the  sound 
wave  does  not  occur.   Thus,  only  phase  weighting  for  geometry 
is  required,  and  only  those  results  will  be  presented.   The 
homogeneous  case  provides  a  needed  baseline  with  which  to 
judge  the  effects  of  model-based  signal  processing  when 
inhomogeneous  case  data  is  studied. 

Since  case  HMG1  includes  the  effects  of  the  ocean  com- 
munication channel,  a  carrier  frequency  is  used  to  convey 
the  modulation  through  the  channel.   In  case  HMG1  and  case 
INHMG1  the  carrier  frequency  was  set  to  5.0  kHz.   The  array 
interelement  spacing  was  adjusted  accordingly  to  maintain 
the  one-half  wavelength  spacing  needed  to  prevent  grating 
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lobes  in  the  array  beam  patterns.   Note  that  the  interelement 
spacing  was  set  using  the  highest  frequency  expected  for 
either  the  CW  pulse  or  the  LFM  pulse  waveforms,  and  was 
kept  constant  for  both  the  HMG1  and  INHMG1  case  studies. 
The  highest  frequency  occurred  for  the  CW  pulse  case,  and  was 
equal  to  6000  Hz.   The  estimates  of  the  time  delay  and 
Doppler  shift  in  subprogram  SGNLGN  were  set  equal  to  the 
actual  values  to  maximize  the  correlator  output. 

The  system  parameters  defining  case  HMG1  are  as  follows: 

-  Array  Parameters 

Number  of  array  elements,  M  =  N  =  5 

Array  Element  Spacing,  d   =  ^   =  0.1229  meters 

-  Medium  Parameters 

Speed  of  Sound,  c   =  1475  meters/second 

Actual  line  of  sight 

time  delay  t   =  2.033898  seconds 

Actual  Doppler  shift  <PA   =  0.0  Hertz 

-  System  Geometry  (See  Figure  1) 

Depth  of  Transmit  Array  y   =  1000.0  meters 

Depth  of  Receive  Array  y   =  2500.0  meters 

Cross  Range  x  -x   =  500.0  meters 

r  r    o 

Line  of  sight  range  |r-r  |  =  3000.0  meters 

Eight  plots  were  generated  for  case  HMG1.   The  first  four 
plots,  Figures  29  through  32,  were  produced  using  the 
rectangular  envelope  CW  pulse  transmit  signal.   The  remain- 
ing four  plots,  Figures  33  through  36  were  produced  by 
the  rectangular  envelope  LFM  pulse  transmit  signal.   Figures 
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29,  31,  33  and  35  show  receiver  performance  when  no  phase 
weighting  is  done.   These  plots  indicate  that  performance 
will  be  very  poor  when  the  output  signals  from  the  array 
elements  are  not  cophased.   Again,  it  is  significant  that 
the  Pd  in  this  case  is  approximately  equal  to  the  value  of 
the  Pfa  parameter  used  in  obtaining  the  plot.   That  is, 
the  total  array  output  signal  must  be  very  small  because 
of  destructive  interference  of  the  array  element  signals, 
and  the  only  detections  appear  to  be  due  to  the  noise. 

Figures  30,  32,  34  and  36  show  performance  when  phase 
weighting  is  done  to  compensate  for  system  geometry.   These 
plots  indicate  receiver  performance  can  be  made  to  approxi- 
mate the  theoretical  predicted  performance  if  the  individual 
array  element  signals  can  be  cophased  to  provi.de  constructive 
interference  at  the  output  of  the  array  processor. 

Figures  29,  30,  33  and  34  show  receiver  performance  at 
a  Pfa  of  0.1.   Figures  31,  32,  35  and  36  show  the  effect  of 
decreasing  the  Pfa  to  0.01.   The  effect  of  decreasing  the 
Pfa  on  the  simulated  performance  follows  the  theoretically 
predicted  performance  and  provides  further  confirmation  that 
the  computer  implementation  of  the  model  is  accurate. 

The  important  conclusion  from  case  HMG1  is  dramatic 
improvement  in  receiver  performance  is  possible  provided 
that  enough  information  about  system  geometry  exists  to 
generate  phase  weights  that  will  properly  cophase  the  array 
element  output  signals. 
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D.   CASE  INHMG1 

The  input  signal  data  for  case  INHMG1  is  produced  in 
the  same  manner  as  case  HMG1.   The  transmit  signal  waveforms, 
system  geometry,  and  transmit  and  receive  arrays  are  identi- 
cal.  The  only  difference  between  the  two  cases  is  the 
model  of  the  ocean  medium.   For  case  INHMG1 ,  the  ocean  is 
considered  to  be  an  inhomogeneous  medium  with  respect  to  the 
speed  of  sound  propagation  through  the  water.   The  speed  of 
sound  is  assumed  to  vary  linearly  with  depth  (or  y-coordinate) 
only.   That  is,  the  speed  of  sound  is  not  identical  at  the 
transmit  and  receive  arrays.   In  this  space-variant  model 
of  the  ocean  medium,  refraction  or  ray  bending  of  the  sound 
wave  does  occur,  and  the  model  of  the  ocean  medium  due  to 
Ziomek  [Ref.  2]  predicts  what  the  phase  shifts  in  array  ele- 
ment output  electrical  signals  will  be.   Using  this  knowledge 
of  the  physics  of  acoustic  wave  propagation,  phase  weights 
are  generated  that  attempt  to  compensate  for  the  refraction 
of  the  acoustic  wave  along  the  ray  path,  and  the  system 
geometry.   To  determine  the  effect  of  this  model-based  signal 
processing  algorithm,  the  inhomogeneous  data  was  also 
processed  using  phase  weights  that  compensate  just  for  sys- 
tem geometry.   The  difference  in  simulated  receiver  perfor- 
mance between  these  two  situations  will  graphically  show  the 
effectiveness  of  the  model-based  signal  processing  approach. 

The  system  parameters  defining  case  INHMG1  are  as 
follows : 
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-  Array  Parameters 

Number  of  array  elements,     M  =  N  =  5 

Array  Element  Spacinq,        d   =  d     0.1229  meters 
2  r  v  x      y 

-  Medium  Parameters 

Speed  of  Sound,  c   =  1475  meters/second 

Gradient  g   =  .017 

Actual  line  of  sight 

time  delay  t,  =  2.033898  seconds 

Actual  Doppler  shift  cf>A  =  0.0  Hertz 

-  System  Geometry  (See  Figure  1) 

Depth  of  Transmit  Array       y   =  1000.0  meters 

Depth  of  Receiver  Array       y   =  2500.0  meters 

Cross  Range  x  -x   =  500.0  meters 

3  r    o 

Line  of  sight  range      |r  -r  |  =  3000.0  meters 
Twelve  plots  were  generated  for  case  INHMG1.   The  first 
six  plots,  Figures  37  through  42  show  receiver  performance 
for  a  rectangular  envelope  CW  pulse.   The  remaining  six 
plots,  Figures  43  through  48,  are  the  result  of  a  rectangu- 
lar envelope  LFM  pulse  transmit  signal.   Differences  in 
receiver  performance  due  to  the  form  of  the  transmit  signal 
are  measurable,  but  not  significant. 

Figures  37,  40,  43  and  46  were  generated  when  no  phase 
weighting  is  applied  to  the  array  output  electrical  signals 
Again,  these  figures  indicate  that  the  array  element  output 
electrical  signals  must  be  cophased  to  provide  any  useful 
detection  capability  in  the  receiver. 
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Figures  38,  41,  44  and  47  show  the  receiver  performance 
when  only  the  system  geometry  is  compensated  for  by  the 
phase  weights  used  in  the  array  processor.   These  figures 
indicate  that  in  the  presence  of  the  inhomogeneous  ocean 
medium,  phase  weighting  to  compensate  for  system  geometry 
will  not  achieve  the  receiver  performance  predicted  by 
theory.   That  is,  traditional  beam  steering  is  not  suffi- 
cient for  maximum  receiver  performance,  and  a  large  margin 
for  improved  performance  exists  when  the  effects  of  the 
ocean  medium  on  the  received  signal  can  be  predicted.   Note 
also  that  the  degradation  in  performance  is  consistent  for 
the  different  transmit  waveforms  and  values  of  Pfa. 

Figures  39,  42,  45  and  48  show  the  effect  on  receiver 
performance  when  phase  weights  are  computed  based  on  both 
the  system  geometry  and  the  refraction  of  the  acoustic  wave 
in  the  ocean  medium.   The  state  of  the  logical  variable 
DMEDIA  indicates  if  the  phase  weights  contain  the  correction 
for  deterministic,  inhomogeneous  medium  wave  propagation 
effects.   As  expected,  the  simulated  receiver  performance 
conforms  to  the  performance  predicted  by  theory,  and  is 
nearly  identical  to  that  of  case  TEST  when  all  medium 
effects  have  been  corrected  for. 

However,  the  really  significant  result  is  seen  by  com- 
paring Figures  38,  41,  44  and  47  with  Figures  39,  42,  45  and 
48.   The  only  difference  between  these  plots  is  in  use  of 
phase  weights  that  compensate  for  the  phase  shifts  due  to 
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propagation  through  an  inhomogeneous  ocean  medium,  and  the 
model  of  the  inhomogeneous  ocean  medium  is  based  on  the 
physics  of  acoustic  wave  propagation.   If  one  assumes  a 
marginal  detection  probability  of  50%  (Pd  =  0.5),  the 
improvement  in  receiver  performance  due  to  the  application 
of  a  model-based  signal  processing  approach  for  this  particu- 
lar simulation  geometry  and  ocean  medium  characteristic 
varies  from  14  dB  to  18  dB  depending  on  the  type  of 
transmit  waveform  or  value  of  Pfa.   The  conclusion  is  that 
the  physics  of  wave  propagation  can  be  used  to  develop 
signal  processing  algorithms  having  significant  impact  on 
the  performance  of  a  receiver  designed  to  process  complex 
envelope,  electrical  signals  from  an  array  of  point  source 
elements . 
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V.   CONCLUSIONS  AND  RECOMMENDATIONS 

The  original  objectives  of  this  research  have  been  accom- 
plished.  The  computer  simulation  of  a  correlator/matched 
filter  receiver  to  process  data  from  an  array  of  electro- 
acoustic  transducers  has  been  developed  and  validated.   The 
concept  of  model-based  signal  processing  was  applied  to  the 
development  of  an  array  signal  processing  algorithm,  and  was 
shown  to  have  a  marked  impact  on  the  capability  of  the 
receiver  to  detect  the  presence  of  a  signal  in  zero-mean, 
additive,  white  Gaussian  noise. 

The  computer  simulation  of  the  receiver  was  validated  by 
a  direct  comparison  between  the  measured  performance  of  the 
simulation  and  the  performance  predicted  by  theory  when  all 
array  element  output  electrical  signals  are  precisely 
cophased.   The  agreement  between  predicted  and  measured 
performance  was  close  but  not  exact  since  the  noise  assump- 
tions made  in  developing  the  theoretical  performance  could 
not  be  precisely  duplicated  by  the  computer  simulation. 
However,  the  close  agreement  with  theoretical  performance 
was  considered  to  be  a  validation  of  the  receiver  simulation 
computer  program. 

Test  cases  were  designed  and  run  that  permitted  relative 
changes  in  receiver  performance  to  be  measured.   In  this 
way,  the  effectiveness  of  the  model-based  signal  processing 
approach  could  be  quantified.   For  the  simulation  parameters 
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and  transmit  signals  used  in  this  study,  compensating  for 
the  refraction  of  the  acoustic  field  in  the  ocean  medium 
was  found  to  improve  receiver  performance  over  traditional 
beam  steering  by  at  least  14  dB  when  measured  at  the  point 
of  50%  probability  of  detection. 

Suggestions  for  further  study  of  the  receiver  model 
include : 

-  Obtain,  or  develop,  a  pseudorandom  noise  generator 
that  more  closely  approximates  the  noise  model 
assumptions  implicit  in  the  derivation  of  the 
theoretical  receiver  performance,  and  revalidate 
the  computer  simulation. 

-  Simulate  more  realistic  noise  processes,  such  as 
colored  noise,  and  determine  the  impact  on  receiver 
performance. 

-  Try  decision  rules  other  than  Neyman-Pearson,  such 

as  the  minimum  average  probability  of  error  criterion 
to  measure  probability  of  error  performance  for  an 
underwater,  acoustic  data  communication  system. 

-  Experiment  with  other  functional  forms  for  the 
transmit  signal,  and  measure  the  effect  on  perfor- 
mance when  errors  in  the  estimate  of  time  delay  and 
Doppler  shift  exist. 
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